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Abstract 

This document describes a Fortran 95 package for carrying out DGLAP evolu- 
tion and other common manipulations of parton distribution functions (PDFs). The 
PDFs are represented on a grid in x-space so as to avoid limitations on the func- 
tional form of input distributions. Good speed and accuracy are obtained through 
the representation of splitting functions in terms of their convolution with a set of 
piecewise polynomial basis functions, and Runge-Kutta techniques are used for the 
evolution in Q. Unpolarised evolution is provided to NNLO, including heavy-quark 
thresholds in the MS scheme, and longitudinally polarised evolution to NLO. The 
code is structured so as to provide simple access to the objects representing splitting 
functions and PDFs, making it possible for a user to extend the facilities already 
provided. A streamlined interface is also available, facilitating use of the evolution 
part of the code from F77 and C/C++. 
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Program Summary 



Title of program: HOPPET 
Version: 1.1 
Catalogue identifier: 

Program obtainable from: http://projects.hepforge.org/hoppet/ 



Distribution format: compressed tar file 

E-mail: salam@lpthe.jussieu.fr, rojo01pthe.jussieu.fr 

License: GNU Public License 
Computers: all 
Operating systems: all 
Program language: Fortran 95 
Memory required to execute: < 10 MB 
Other programs called: none 
External files needed: none 

Number of bytes in distributed program, including test data etc.: ~ 270 kB 

Keywords: unpolarised and longitudinally polarised parton space- like distribution func- 
tions (PDFs), DGLAP evolution equations, x-space solutions. 

Nature of the physical problem: Solution of the DGLAP evolution equations up to NNLO 
(NLO) for unpolarised (longitudinally polarised) PDFs, and provision of tools to facilitate 
manipulation (convolutions, etc.) of PDFs with user-defined coefficient and splitting func- 
tions. 

Method of solution: representation of PDFs on a grid in x, adaptive integration of split- 
ting functions to reduce them to a discretised form, obtaining fast convolutions that are 
equivalent to integration with an interpolated form of the PDFs; Runge Kutta solution of 
the Q evolution, and its caching so as to speed up repeated evolution with different initial 
conditions. 

Restrictions on complexity of the problem: PDFs should be smooth on the scale of the 
discretisation in x. 

Typical running time: a few seconds for initialisation, then ~ 10 ms for creating a tabula- 
tion with a relative accuracy of 10~^ from a new initial condition (on a 3.4 GHz Pentium IV 
processor). Further details in Sect. 
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1 Introduction 

There has been considerable discussion over the past years (e.g. [Il[2],[3llll[5l[6l[7l[8], 
[HI [in]) of numerical solutions of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) 
equation [TT] for the Quantum Chromodynamics (QCD) evolution of parton distribution 
functions (PDFs). 

The DGLAP equation [11] is a renormalisation group equation for the quantity qi{x, Q"^), 
the density of partons of type (or flavour) i carrying a fraction x of the longitudinal mo- 
mentum of a hadron, when resolved at a scale Q. It is one of the fundamental equations 
of perturbative QCD, being central to all theoretical predictions for hadron-hadron and 
lepton-hadron colliders. 

Technically, it is a matrix integro-differential equation, 

dqi{x,Q^) asiQ^) f'dz 



^ ^P.,{z,a,{Q'))q,{^^,Q') , (1) 



whose kernel elements Pij{z,Q ) are known as splitting functions, since they describe the 
splitting of a parton of kind j into a parton of kind i carrying a fraction z of the longitudinal 
momentum of j. The parton densities themselves qi{x, Q"^) are essentially non-perturbative, 
since they depend on physics at hadronic mass scales < 1 GeV, where the QCD coupling 
is large. On the other hand the splitting functions are given by a perturbative expansion 
in the QCD coupling ^^(Q^). Thus given partial experimental information on the parton 
densitie^ — for example over a limited range of Q, or for only a subset of parton flavours 



^Of course it is not the parton densities, but rather structure functions, which can be derived from 
them perturbatively, that are measured experimentally. 
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— the DGLAP equations can be used to reconstruct the parton densities over the full 
range of Q and for all flavours. 

The pivotal role played by the DGLAP equation has motivated a considerable body of 
literature discussing its numerical solution PjElEllllElinilTllHlE], [in]. There exist two 
main classes of approaches: those that solve the equation directly in x-space and those 
that solve it for Mellin transforms of the parton densities, defined as 



and subsequently invert the transform back to x-space. Recently, a novel approach has been 
proposed which combines advantages of the A^— space and x— space methods [8]. A^— space 
based methods are of interest because the Mellin transform converts the convolution of 
eq. ([1]) into a multiplication, resulting in a continuum of independent matrix differential 
(rather than integro-differential) equations, one for each value of A^, making the evolution 
more efficient numerically. 

The drawback of the Mellin method is that one needs to know the Mellin transforms 
of both the splitting functions and the initial conditions. There can also be subtleties 
associated with the inverse Mellin transform. The x-space method is in contrast more 
flexible, since the inputs are only required in x-space; however it is generally considered to 
less efficient numerically, because of the need to carry out the convolution in eq. ([I]). 

To understand the question of efficiencies one should analyse the number of operations 
needed to carry out the evolution. Assuming that one needs to establish the results of 
the evolution at A'^, values of x, and Aq values of Q, one essentially needs O {N^Nq) 
operations with an x-space method, where the factor comes from the convolutions. In 
the Mellin-space method, one needs O (N^NqM) operations, where M is the number of 
points used for Mellin inversion. One source of drawback of the x-space method is that, 
nearly always, A^^^ ~ Inl/xmm and so the method scales as the square of Inl/xmm, where 
the Mellin method is linear (and M can be kept roughly independent of Xmin) • 

The other issue relates to how one goes to higher numerical integration and interpolation 
orders. In x-space methods one tends to choose x values that are uniformly distributed 
(be it in In 1/x or some other more complex function) — this limits one to higher-order 
extensions of the Trapezium and Simpson-rule type integrations, whose order in general 
is np — 1 where Up is the number of points used for the integration. The precision of 
the integration is given by (5x)"^ where 6x is the grid spacing. Higher Up improves the 
accuracy, but typically Up can not be taken too large because of large cancellations between 
weights that arise for large Up. In the Mellin method one is free to position the M points 
as one likes, and one can then use Gaussian type integration O [3], [9]; using Hp points 
one manages to get a numerical order 2np — 1, i.e. accuracy {6x)'^^'' , and furthermore the 
integration weights do not suffer from cancellations at large Up, allowing one to increase 
Hp, and thus the accuracy, quite considerably. 

Despite it being more difficult to obtain high accuracy with x-space methods, their 
greater flexibility means that they are widespread, serving as the basis of the well-known 
QCDNUM program P, and used also by the CTEQ [12] and MRST/MRSW [I3] global 




(2) 
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fitting collaborations. Higher-order methods in x-space have been developed in [21 HJ [6l [71 
[H], and more recently have been incorporated also in QCDNUM. 

HOPPET, the program presented here, uses higher-order methods both for the x-integrations 
and Q evolution. It combines this with multiple grids in x-space: a high-density grid at 
large x where it is hardest to obtain good accuracy, and coarser grids at smaller x where 
the smoothness of the PDFs facilitates the integrations. One of the other crucial features 
of the program is that it pre-calculates as much information as possible, so as to reduce 
the evolution of a new PDF initial condition to a modest set of addition and multiplication 
operations. Additionally, the program provides access to a range of low and medium-level 
operations on PDFs which should allow a user to extend the facilities already provided. 

The functionality described in this article has been present in hoppet's predecessors 
for several years (they were available on request), but had never been documented. Those 
predecessors have been used in a number of different contexts, like resummation of event 
shapes in DIS [13], automated resummation of event shapes [TS], studies of resummation in 
the small-x limit [TB], and in a posteriori inclusion of PDFs in NLO final-state calculations 
[T7t [18], as well as used for benchmark comparisons with Pegasus [3] in [19]. Since the 
code had not hitherto been released in a documented form, it is the authors' hope that 
availability of this documentation may make the package somewhat more useful. 

This manual is structured as follows: Sect. [2] briefly summarises the perturbative QCD 
ingredients contained in hoppet, while Sect. [3] describes the numerical techniques used 
to solve the DGLAP equation. Sects. [3HZI present in detail the capabilities of the hoppet 
package with its general F95 interface, with emphasis on those aspects that can be adapted 
by a user to tailor it to their own needs. Sect. [8] describes a streamlined interface to 
HOPPET which embodies its essential capabilities in a simple interface available a variety of 
programming languages: F77 and C/C-I--I-. Finally, Sect. [3 presents a detailed quantitative 
study of the performance of hoppet, and in the final section we conclude. A set of 
appendices contain various example programs, both for the general and the streamlined 
interfaces, a reference guide with the most important hoppet modules, details on technical 
aspects and a set of useful tips on Fortran 95. 

A reader whose interest is to use hoppet to perform fast and efficient evolution of 
PDFs may wish to skip Sects. [H to [7] and move directly to Sect. [HI which describes the 
user-friendlier streamlined interface, accessible from F95, F77 and C/C-I--I-, and which 
contains the essential functionalities of hoppet. He/she is also encouraged to go through 
the various example programs which contain detailed descriptions and explanations. On 
the other hand, a reader interested in the more flexible and general functionalities of 
hoppet, should also consult Sects. [3|to[71 

Note that throughout this documentation, a PDF refers always to a momentum density 
rather than a parton density, that is, when we refer to a gluon, we mean xg{x) rather than 
g{x), the same convention as used in the LHAPDF PDF library [20] . 
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2 Perturbative evolution in QCD 



First of all we set up the notation and conventions that are used through hoppet. The 
DGLAP equation for a non-singlet parton distribution reads 

(3) 

Note that the related variable t = InQ^ is also used through hoppet. The splitting 
functions in eq. ([3]) are known up to NNLO in the unpolarised case pTl [22] : 

P(., (Q^) ) = P(°) (.) + ^iMpd) (,) + (^^iM^ ' p(2) (,) , (4) 

and up to NLO in the polarised case. The generalisation to the singlet case is straightfor- 
ward, as it is the generalisation of eq. ([3]) to the case of time-like evolutioiJl, relevant for 
example for fragmentation function analysis, where partial NNLO results are also available 

m- 

As with the splitting functions, all perturbative quantities in hoppet are defined to 
be a coefficient of agjl-n. The one exception is the /5-function coefficients of the running 
coupling equation: 

= /? («. (Q')) = -as(/?oa. + Aa' + /52«D • (5) 

The evolution of the strong coupling and the parton distributions can be performed 
in both the fixed flavour-number scheme (FFNS) and the variable flavour-number scheme 
(VFNS). In the VFNS case we need the matching conditions between the effective theories 
with rij and rij + 1 light flavours for both the strong coupling (Q^) and the parton 
distributions at the heavy quark mass threshold m^. 

These matching conditions for the parton distributions receive non-trivial contributions 
beginning at NNLO. For light quarks qi^i of flavour i (quarks that are considered massless 
below the heavy quark mass threshold m^) the matching between their values in the 
and 77, J -(- 1 effective theories reads: 

where i = 1, . . . n/, while for the gluon distribution and the heavy quark PDF one has 

^ The general structure of the relation between space-like and time-like evolution and splitting functions 
has been investigated in [23 l [24 l [25 l [26 l [27] . 
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a coupled matching condition: 
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with qh = qhi and the singlet PDF S(a;, Q^) is defined in Table [TJ The NNLO matching 
coefficients were computed in |28j^ . Notice that the above conditions will lead to small 
discontinuities of the PDFs in its evolution in Q^, which are cancelled by similar matching 
terms in the coefficient functions resulting in continuous physical observables. In particular, 
the heavy quark PDFs start from a non-zero value at threshold at NNLO, which sometimes 
can even be negative. 

The corresponding NNLO relation for the matching of the coupling constant at the 
heavy quark threshold m\ is given by 

ai"'«\m?J=ai»''K) + cJ^^^^y . (8) 

where the matching coefficient C2 was computed in |29j . 

The default basis for the PDFs, called the human representation in hoppet, is such 
that the entries in an array pdf (-6:6) of PDFs correspond to: 

t = — 6 , 6 = — 5 , c = — 4 , s = — 3 , u = — 2 , J = — 1 , 

= 0, (9) 
(i = l,M = 2,s = 3,c = 4,6 = 5,t = 6. 



This representation is the same as that used through the LHAPDF library [20]. However, 
this representation leads to a complicated form of the evolution equations. The splitting 
matrix can be simplified considerably (made diagonal except for a 2 x 2 singlet block) by 
switching to a different ffavour representation, which is named the evln representation, for 
the PDF set, as explained in detail in [301 El]. This representation is described in Table HJ 

In the evln basis, the gluon evolves coupled to the singlet PDF S, and all non-singlet 
PDFs evolve independently. Notice that the representations of the PDFs are preserved 
under linear operations, so in particular they are preserved under DGLAP evolution. The 
conversion from the human to the evln representations of PDFs requires that the number 
of active quark ffavours Uf is specified by the user, as described in Sect. I5.1.2[ 

In HOPPET unpolarised DGLAP evolution is available up to NNLO in the MS scheme, 
while for the DIS scheme only evolution up to NLO is available. For polarised evolution 
only the MS scheme is available. The variable f actscheme takes different values for each 
factorisation scheme: 



^The authors are thanked for the code corresponding to the calculation. 
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name 




-Q...-(nf + l) 




Qi 


-nf...-2 


5'NS,j 


(Qi - Qi) - (Qi - Qi) 


-1 




Y.7=i{Qj-Qj) 





g 


gluon 


1 


S 


Y.7=i{Qj + Qj) 


2. . .n/ 




{Qi + Qi) - {Qi + Qi) 


{nf + l)...Q 




Qi 



Table 1: 

The evolution representation (called evln in hoppet) of PDFs with rif active quark 
flavours in terms of the human representation. 



fact scheme 


Evolution 


1 


unpolarised MS scheme 


2 


unpolarised DIS scheme 


3 


polarised MS scheme 



Note that mass thresholds are currently missing in the DIS scheme. 

3 Numerical techniques 

We briefly introduce now the numerical techniques used to perform parton evolution: the 
discretisation of PDFs and their convolutions with splitting functions on a grid in x, and 
the subsequent DGLAP evolution in . 

The first aspect that we discuss is how to represent PDFs and associated convolutions in 
terms of an interpolating grid in x— space. Note that this technique can be applied to any 
quantity that appears through convolutions, so it is not restricted to parton distributions. 
Then we will describe how to obtain the solution of the DGLAP evolution equations. 

3.1 Higher order matrix representation 

Given a set of grid points Ua = IuI/xq, labelled by an index a and (for later con- 
venience) a uniform grid spacing, Ua = aSy, one can approximate a parton distribution 
function xq (x, t) by interpolating the PDF at the grid pointsQ 

xq{y = \nl/x,t) = ^Wa{y)qa{t) , (10) 



'^From the numerical point of view, it is advantageous to interpolate xq rather than q itself because the 
is in general smoother. 
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where Waiy) are the interpolation weights, we have defined 

qciit) = XaqiVa^t) , (11) 

and the sum over a runs over n + 1 points in the vicinity of y for an interpolation order 
n. Note that Greek indices represent the y dimension, while Roman indices are used to 
represent the fiavour dimension. 

The convolution of a single-fiavour PDF with a splitting function P{z, t) can be written 

as 

{P^q){y,t) = J2^a{y){P®q)a{t) , (12) 

a 

where we have replaced the convolution by its grid representation, 

{P^qUt) = J2P-p{t)qp{t), (13) 

13 

where (3 runs over O {N^) points of the grid and we have defined 

Pa/sit) = f dzP{z,t)wp{y^ + \nz) . (14) 

Note that the piecewise interpolating polynomials that we use are not smooth at the 
grid points. The reason for this is that this type of interpolation strategy is much more 
suited for the numerical integration required in convolutions than other strategies which 
impose continuity on the interpolating polynomials, like for example splines, which are in 
general slower and less accurate for the same number of grid points. 

A virtue of having a uniform grid in y = In 1/x is that the interpolation functions 
can be arranged to have a structure wpiy) = wp{y — yp), so that Pap just depends on 
a — (3, and can be rewritten Va-p- A slight subtlety arises at large x, where if one writes 
'^piy) = '^phj ~ Up) is effectively assuming an interpolation that uses identically zero 
interpolation points for x > 1 [6l [10], even though xq{x) is not formally defined for x > 1. 
In practice this is often not too important (because PDFs drop rapidly towards x = 1), 
but we shall include two options: one that uses effective zero-points beyond x = 1 and 
one that ensures that the interpolation is based only on the physically valid domain of the 
PDFs. 

These two choices are represented in Fig. [H For an interpolation of order n (that 
is, which uses information from n + 1 grid points), the option of using only points with 
X < 1 is denoted by order = n. This has the consequence that for jS < n we cannot write 
Pai3 = 'Pa- 13, and so must explicitly store O {N^n) distinct Pa/3 entries. The option of using 
artificial (zero-valued) points at x > 1 is denoted by order = — n, and does allow us to 
write Pa/3 = Va-zs (thus we store only O (N^) entries), with Va~/3 = for a < /3. 
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X<1 X=l X>1 

ORDER = 5 : 



ORDER = 5 



a: 987654321 = 



PDF = 



ORDER = -5 



ORDER = -5 



Figure 1: The different strategies (+ve and — ve order) for interpolation of the grid near 
X = 1. The dotted (blue, red) boxes indicate two regions in which we illustrate the 
interpolation of the PDF, while the (blue, red) lines with arrows indicate the corresponding 
range of grid points on which the interpolation is based. 

3.2 Evolution operators 

The DGLAP evolution equation, eq. ([3]), is easily approximated in terms of its grid repre- 
sentation by 

'-^ = ^^PMt), (15) 

where for a general value of a the sum over f] extends over O {a + | order |) points of 
the grid. Introducing Mapito) = Sajs as initial condition at some initial scale to, one can 
alternatively solve 

— -qI = 2^ Po,^[t)M^p[t) . (16) 

7 

Then the evolved parton distribution at the grid points is given by 

q^{t) = Y,Mc.p{t)qp{to). (17) 

We refer to Ma/sit) as the evolution operator. From a practical point of view, we will solve 
eqs. f fTSjl and ([T6|l with higher order iterative Runge-Kutta methods, as described in Sect. 

Em 

A further simplification occurs if one can rewrite the splitting functions P as a trans- 
lation invariant object. Using the properties discussed in the previous subsection, we can 
rewrite P^/? = Va-fs, and then similarly one can rewrite Map = Aia-i3, and it is as simple 
to determine Maf}{t) as it is to determine the evolution of a single vector q^, i.e. one just 
evolves a single column, /9 = 0, of M^pit). 
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4 Single- flavour grids and convolutions 



HOPPET is written in Fortran 95 (F95). This has considerable advantages compared to 
F77, as will be seen in the discussion of the program, though it does lack a number of fully 
object-oriented features and this sometimes restricts the scope for expressiveness. For- 
tran 95 perhaps not being the best known language in the high-energy physics community, 
occasionally some indications will be give to help the reader with less-known language 
constructs, with further information in Appendix [El 

All routines described in this section need access to the convolution module, which 
can either be obtained directly by adding a 

use convolution 

statement at the beginning of the relevant subprogram (before any implicit none or 
similar declarations). Alternatively, as with the rest of the routines described in this 
documentation, it can be accessed indirectly through the hoppet_vl module 

use hoppet_vl 

Unless you are delving into the innards of hoppet, the latter is more convenient since it 
provides access to everything you are likely to need. 

4.1 Grid definitions (grid_def) 

The grid (in y) is the central element of the PDF evolution. Information concerning the 
grid is stored in a derived type grid_def : 

type(grid_def ) : : grid 

call InitGridDef (grid , dy=0 . l_dp , ymax=10 . 0_dp , order=5) 

This initialises a grid between x = 1 and down to a: = e~^^^, with uniform grid spacing in 
?/ = In 1/x of dy=0. 1, with a grid which uses order 5 interpolation with only x < 1 points. 
The user can modify this choice to better suit his/her particular needs, as explained in 
Sect. 13. 1[ One notes the use of keyword arguments — the keywords are not mandatory in 
this case, but have been included to improve the legibility. Having defined a grid, the user 
need not worry about the details of the grid representation. 

In line with the convention set out in the Fortran 90 edition of Numerical Recipes [32] 
we shall use _dp to indicate that numbers are in double precision, and real (dp) to declare 
double precision variables. The integer parameter dp is defined in the module types (and 
available indirectly through module hoppet_vl). 

It is often useful to have multiple grids, with coarser coverage at small x and finer 
coverage at high x, to improve the precision the convolutioii^l at large-x without introducing 
an unnecessarily large density of points at small-x. To support this option, we can first 
define an array of sub-grids, and then use them to initialise a combined grid as follows: 

The reason that denser grids are required at largc-x is that if a typical parton distributions goes as 

hmg(a;,Q2) ^ ^ (18) 
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type (gr id_d.ef ) 



grid, subgridsO) 



! define the various sub-grids 

call InitGridDef (subgrids(l) ,dy=0.30_dp, ymax=10.0_dp, order=5) 
call InitGridDef (subgrids (2) ,dy=0. 10_dp, ymax= 2.0_dp, order=5) 
call InitGridDef (subgrids(3) ,dy=0.03333_dp, ymax= 0.6_dp, order=5) 
! Smaller dy at small ymax / large xmin 

! put them together into a single combined grid 
call InitGridDef (grid, subgrids, locked= . true . ) 

When combining them, the locked=.true. option has been specified, which ensures that 
after any convolution, information from the finer grids is propagated into the coarser ones. 
This places some requirements on the grid spacings, notably that a coarse grid have a 
spacing that is a multiple of that of the next finest grid. If the requirements are not 
satisfied by the subgrids that have been provided, then new similar, but more suitable 
subgrids are automatically generated. 

Note that the two kinds of invocation of InitGridDef actually correspond to different 
(overloaded) subroutines. The Fortran 95 compiler automatically selects the correct one 
on the basis of the types of arguments passed. 

Though only grids that are uniform in y have been implemented (and the option of a 
simultaneous combination of them), nearly all of the description that follows and all code 
outside the convolution module are independent of this detail, the only exception being 
certain statements about timings. Therefore were there to be a strong motivation for an 
alternative, non- uniform grid, it would suffice to modify the convolution module, while 
the rest of the library (and its interfaces) would remain unchanged. 



4.2 £c-space functions 

Normal x-space functions (such as PDFs) are held in double precision arrays, which are to 
be allocated as follows 

real(dp), pointer :: xgluon(:) 
call AllocGridQuant(grid,xgluon) 

Note that for this to work, xgluon(:) should be a pointer, and not just have the 
allocatable attribute. To deallocate a grid quantity, one may safely use the F95 deallocate 
command. 

Since xgluon( : ) is just an array, it carries no information about the grid. Therefore 
to set and access its value, one must always provide the information about the grid. This 
is not entirely satisfactory, and is one of the drawbacks of the use of F95. 

then its logarithmic derivative with respect to x is divergent, 

dlnq{x,Q'^) -mx 
lim — = lim *■ — oo , (19) 

x^i omx x^i 1 — x 

and therefore to maintain the relative accuracy of the evolution, grids with denser coverage at large- cc are 
required. 
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There are a number of ways of setting a grid quantity. Suppose for example that we 
have a function 

function exajiiple_xgluon(y) 

use types ! ! defines "dp" (double precision) kind 

implicit none 

real (dp), intent (in) :: y 
real (dp) : : x 
X = exp(-y) 

excmiple_xgluon = 1.7_dp * x**(-0.1_dp) * (l-x)**5 !! returns xg(x) 
end function excmiple_xgluon 

which returns the gluon momentum density xg{x) (cf. Sect. 13. ip . Then we can call 

call InitGr idQuant (gr id , xgluon , exainple_xgluon) 

to initialise xgluon with a representation of the return value from the function example_xgluon. 
Alternative methods for initialising grid quantities are described in Appendix O 
To then access the gluon at a given value of |/ = In 1/x, one proceeds as follows 

real (dp) : : y, xgluon_at_y 
y = 5.0_dp 

xgluon_at_y = EvalGr idQuant (grid, xgluon, y) ! ! again this returns xg(x) 

Note that again we have to supply the grid argument to EvalGridQuant because the 
xgluon array itself carries no information about the grid (other than its size). 

A less efficient, but perhaps more 'object-oriented' way of accessing the gluon is via the 
notation 

xgluon_at_y = xgluon .aty. (y. with. grid) 

There also exists an .atx. operator for evaluating the PDF at a given x value. Many 
of these procedures and operators are overloaded so as to work with higher- dimensional 
arrays of grid quantities, for example a multi-flavour PDF array pdf (:,:). The flrst index 
will always correspond to the representation on the grid, while the second index would 
here indicate the flavour. 

Note that arithmetic operators all have higher precedence than library-deflned operators 
such as . aty . ; accordingly some ways of writing things are more efficient than others: 

xgluon_at_y_times_2 = 2 * xgluon .aty. (y. with. grid) ! very inefficient 
xgluon_at_y_times_2 = 2 * (xgluon .aty. (y. with. grid)) ! fairly efficient 
xgluon_at_y_times_2 = 2 * EvalGridQuant(grid, xgluon, y) ! most efficient 

In the flrst case the whole of the array xgluon is multiplied by 2, and then the result is 
evaluated at y, whereas in the second and third cases only the result of the gluon at y is 
multiplied by 2. 

4.3 Grid convolution operators 

While it is relatively straightforward internally to represent a grid-quantity (e.g. a PDF) as 
an array, for convolution operators it is generally useful to have certain extra information. 
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Accordingly a derived type has been defined to hold a convolution operator, and routines 
are provided for allocation and initialisation of splitting functions. The following example 
describes how the gg LO splitting function would be used to initialise the corresponding 
convolution operator: 

type(grid_conv) : : xPgg 

call AllocGridConvCgrid.xPgg) 

call InitGridConvCgrid.xPgg, xPgg_func) 

where the Pgg splitting function is provided in the form of the function xPgg_func. Note 
that this function must return xPgg{x): 

! returns various components of exp(-y) P_gg (exp(-y)) 
real(dp) function xPgg_func(y) 
use types 

use convolution_communicator ! provides cc_piece, and cc_REAL, . . . 

use qcd ! provides CA, TR, nf , ... 

implicit none 

real (dp), intent (in) :: y 

real (dp) : : x 

X = exp(-y) ; xPgg_func = zero 

if (cc_piece == cc_DELTA) then ! Delta function term 

xPgg_func = (11*CA - 4*nf *TR)/6.0_dp 
else 

if (cc_piece == cc_REAL .or. cc_piece == cc_REALVIRT) & 

& xPgg_func = 2*CA* (x/ (one-x) + (one-x)/x + x*(one-x)) 
if (cc_piece == cc_VIRT .or. cc_piece == cc_REALVIRT) & 

& xPgg_func = xPgg_func - 2*CA*one/ (one-x) 
xPgg_func = xPgg_func * x ! remember to return x * Pgg 

end if 
end function xPgg_func 

To address the issue that convolution operators can involve plus-distributions and delta 
functions, the module convolution_coinmunicator contains a variable cc_piece which 
indicates which part of the splitting function is to be returned — the real, virtual, real + 
virtual, or 5-function pieces. 

The initialisation of a grid_conv object uses an adaptive Gaussian integrator (a variant 
of CERNLIB's dgauss) to calculate the convolution of the splitting function with trial 
weight functions. The default accuracy for these integrations is 10~^. It can be modified 
to value eps with the following subroutine call 

call SetDef aultConvolutionEps(eps) 

which is to be made before creating the grid_def object. Alternatively, an optional eps 
argument can be included in the call to InitGridDef as follows: 

type(grid_def ) : : grid 

real (dp) : : eps 

[ ... set eps ... ] 

call InitGridDef (grid, dy=0 . l_dp , ymax=10 . 0_dp , order=3 , eps) 
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Note that eps is just one of the parameters affecting the final accuracy of convolutions. 
In practice (unless going to extremely high accuracies) the grid spacing and interpolation 
scheme are more critical. 

Having allocated and initialised a xPgg splitting function, we can go on to use it. For 
example: 

real (dp), pointer :: xPgg_x_xgluon( : ) 

call AllocGridQuaiit(grid,xPgg_x_xgluon) ! ! Allocate memory for result of convolution 
_x_xgluon = xPgg .conv. xgluon ! ! Convolution of xPgg with, xgluon 



Since the return value of xPgg .conv. xgluon is just an F95 array, one can also write 
more complex expressions. Supposing we had defined also a xPgq splitting function and a 
singlet quark distribution xquark, as well as as2pi = as/ 271, then to first order in we 
could write the gluon evolution through a step dt in In as 

xgluon = xgluon + (as2pi*dt) * ((xPgg .conv. xgluon) + (xPgq .conv. xquark)) 

Note that like . aty . , . conv . has a low precedence, so the use of brackets is important to 
ensure that the above expressions are sensible. Alternatively, the issues of precedence can 
be addressed by using * (also defined as convolution when it appears between a splitting 
function and a PDF) instead of . conv . : 

xgluon = xgluon + (as2pi*dt) * (xPgg*xgluon + xPgq*xquark) 

Note that, for brevity, from now on we will drop the explicit use of x in front of names 
PDF and convolution operator variables. 



4.3.1 Other operations on grid_conv objects 

It is marginally less transparent to manipulate grid_conv types than PDF distributions, 
but still fairly simple: 



call AllocGridConv(grid,Pab) 

call InitGridConv(grid,Pab) 

call InitGridConvCPab , Pcd [ .factor] ) 

call InitGridConvCgrid , Pab , function) 

call SetToZero(Pab) 

call Multiply (Pab, factor) 

call AddWithCoef f (Pab , Pcd [ , coef f ] ) 

call AddWithCoeff (Pab, function) 



Pab memory allocated 
Pab = (opt. alloc) 

Pab = Pcd [*f actor] (opt. alloc) 
Pab = function (opt. alloc) 

Pab = 

Pab = Pab * factor 

Pab = Pab + Pcd [*coeff] 

Pab = Pab + function 



call SetToConvolution(Pab,Pac,Pcb) ! Pab = Pac.conv.Pcb (opt. alloc) 
call SetToConvolution(P(: , :) ,Pa(: , :) ,Pb(: , :)) ! (opt. alloc) 



! P( 

call SetToCommutator(P( : , : ) ,Pa( : , : ) ,Pb( : , 

! P( 
I 



,:) = matmul(Pa( : , : ) ,Pb( : , : ) ) 
)) ! (opt. alloc) 

,:) = matmul(Pa( : , : ) ,Pb( : , : ) ) 
-matmul(Pb( : , : ) ,Pa( : , : ) ) 



call Delete (Pab) 



! Pab memory freed 
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Routines labelled "(opt . alloc . )" allocate the memory for the grid_conv object if the 
memory has not already been allocated. (If it has already been allocated it is assumed to 
correspond to the same grid as any other grid_conv objects in the same subroutine call). 
Some calls require that one specify the grid definition being used (grid), because otherwise 
there is no way for the subroutine to deduce which grid is being used. 

If repeatedly creating a grid_conv object for temporary use, it is important to remem- 
ber to Delete it afterwards, so as to avoid memory leaks. 

Nearly all the routines are partially overloaded so as to be able to deal with one and two- 
dimensional arrays of grid_conv objects as well. The exceptions are those that initialise 
the grid_conv object from a function (arrays of functions do not exist), as well as the 
convolution routines (for which the extension to arrays might be considered non-obvious) 
and the commutation routine which only has sense for matrices of grid_conv objects. 

4.3.2 Derived grid_conv objects 

Sometimes it can be cumbersome to manipulate the grid_conv objects directly, for example 
when trying to create a grid_conv that represents not a fixed order splitting function, but 
the resummed evolution from one scale to another. For such situations the following 
approach can be used 

real (dp), pointer :: probes (:,:) 
type(grid_coiiv) : : Pqg, Pgq, Presult 
integer : : i 

call GetDerivedProbesCgrid, probes) ! get a set of 'probes' 

do i = 1 , size(probes,dim=2) ! carry out operations on each of the probes 

probes (:,i) = Pqg* (Pgq*probes( : , i) ) - Pgq* (Pqg*probes( : , i) ) 
end do 

call AllocGridConv(grid, Presult) 

call SetDerivedConv(Presult, probes) ! Presult = [Pqg, Pgq] 

Here GetDerivedProbes allocates and sets up an array of probe parton distributions. Since 
a single-flavour parton distribution is a one-dimensional array of real (dp), the array of 
probes is a two-dimensional array of real (dp) , the second dimension corresponding to the 
index of the probe. One then carries out whatever operations one wishes on each of the 
probes. Finally with the call to SetDerivedConv, one can reconstruct a grid_conv object 
that corresponds to the set of operations just carried out 

Some comments about memory allocation: the probes are automatically allocated and 
deallocated; in contrast the call to SetDerivedConv(Presult , probes) knows nothing 
about the grid, so Presult must have been explicitly allocated for a specific grid be- 
forehand. 

A note of caution: when one's grid is made of nested subgrids with the locking option 
set to .true., after a convolution of a grid_def object with a parton distribution, the 
coarser grids for the parton distribution are supplemented with more accurate information 
from the finer grids. When carrying out multiple convolutions, this happens after each 
convolution. There is no way to emulate this with a single grid_def object, and the locking 
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would actually confuse the reconstruction of resulting grid_def object. So when the user 
requests the probes, locking is temporarily turned off globally and then reestablished after 
the derived grid_object has been constructed. Among other things this means that acting 
with a derived grid_object will not be fully equivalent to carrying out the individual 
operations separately. In particular the accuracy may be slightly lower (whatever is lost 
due to the absence of intermediate locking). 

5 Multi-flavour grids and convolutions 

The discussion in the previous section about how to represent functions and associated 
convolutions in a general x— space grid holds for any kind of problem involving convolutions, 
even if the examples were given in the context of DGLAP evolution. In this section we shall 
examine the tools made available specifically to address the DGLAP evolution problem. 

5.1 Full- flavour PDFs and flavour representations 

The routines described in this section are available from the pdf _general and pdf .representation 
modules, or via the hoppet_vl general module. 

Full flavour PDFs sets are just like single flavour PDFs except that they have an extra 
dimension. They are represented by arrays, and if you want hoppet to deal with allocation 
for you, they should be pointer arrays. One can allocate a single PDF (two dimensional 
real (dp) array) or an array of PDFs (three-dimensional real (dp) array) 

real (dp), pointer :: PDF (:,:), PDFarray (:,:,: ) 

call AllocPDF(grid, PDF) ! allocates PDF(0 :, -6 : 7) 

call AllocPDF(grid, PDFarray, 0, 10) ! allocates PDFarray(0 :, -6 : 7 , : 10) 

The first dimension corresponds to the grid in the second dimension corresponds to the 
flavour index. Its lower bound is —6, as one would expect. 

What takes a bit more getting used to is that its upper bound is 7. The reason is 
as follows: the flavour information can be represented in different ways, for example each 
flavour separately, or alternatively as singlet and non-singlet combinations. In practice both 
are used inside the program and it is useful for a PDF distribution to have information 
about the representation, and this is stored in PDF( : ,7)[§. 

5.1.1 Human representation. 

When a PDF is allocated it is automatically labelled as being in the human representation, 
described in Sect. [2l Constants with names like if lv_bbar, if lv_g, if lv_b, are deflned 
in module pdf .representation, to facilitate symbolic access to the different flavours. 

^ In the current release of hoppet, in particular, for PDFs in the human representation one has 
PDF( : ,7)=0, while for PDFs in the evln representation, the information on the active number of flavours 
is encoded as nf = (abs (q(2, 7) )+abs (q(3,7)))/ (abs (q(0 , 7) )+abs (q(l ,7) ) ) , so that it is conserved 
under linear combinations. 

However this feature might be modified in future versions of the program. 
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If you are creating a PDF as an automatic array (one whose bounds are decided not 
by the allocation routine, but on the fly), for example in a function that returns a PDF, 
then you should label it yourself as being in the human representation, either with the 
LabelPdf AsHumanCpdf ) subroutine call, or by setting pdf ( : ,7) to zero: 

module pdf _iiiitial_conditioii 

use hoppet_vl; implicit none 
contains 

function unpolarized_dummy_pdf (xvals) result (pdf) 
real (dp), intent (in) :: xvals(:) 
real (dp) :: pdf (size (xvals) , -6 : 7) 

! clean method for labelling as PDF as being in the human representation 
call LabelPdf AsHuman(pdf) 

! Alternatively, by setting everything to zero 
! (notably pdf(:,7)), the PDF representation 
! is automatically set to be human 
pdf(:,:) = 

! iflv_g is pre-defined integer parameter (=0) for symbolic ref . to gluon 
pdf ( : ,if lv_g) = 1.7_dp * xvals** (-0 . l_dp) * (l-xvals)**5 ! Returns x*g(x) 
[. . . set other flavours here . . .] 
end function unpolarized_dummy_pdf 
end module pdf _initial_condition 

The function has been placed in a module so as to provide an easy way for a calling routine 
to have access to its interface (this is needed for the dimension of xvals to be correctly 
passed) . Writing a function such as that above is probably the easiest way of initialising a 
PDF: 

use hoppet_vl; use pdf _initial_condition; implicit none 
type(grid_def ) : : grid 
real (dp), pointer :: pdf ( : , : ) 
[. . .] 

call AllocPDF(grid,pdf ) 

pdf = unpolarized_dummy_pdf (xValues(grid) ) 

There exist a number of other options, which can be found by browsing through src/pdf _general . f 90. 
Of these a sometimes handy one is 

call AllocPDF(grid,pdf ) 

call lnitPDF_LHAPDF(grid, pdf, LHAsub, Q) 
where LHAsub is the name of a subroutine with the same interface as LHAPDF's evolvePDF 




subroutine LHAsub(x,Q,res) 
use types; implicit none 
real(dp), intent(in) :: x,Q 

real(dp), intent(out) :: res(-6:6) ! on output contains flavours -6:6 at x,Q 
[...] ! Note that it should return momentum densities 

end subroutine LHAsub 
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Note that LHAsub should return momentum densities, as happens with the LHAPDF rou- 
tines [20]. 

Having initiahsed a PDF, to then extract it at a given y value, one can either examine 
a particular flavour using the methods described in Sect. 14.21 

real (dp) : : y, gluon_at_y 

gluon_at_y = pdf ( : , if lv_g) .aty. (y. with. grid) 
! OR 

gluon_at_y = EvalGridQuant (grid, pdf (: ,iflv_g) ,y) 

or one can extract all flavours simultaneously 

real(dp) :: pdf _at_y(-6 : 6) 

pdf_at_y = pdf(:,-6:6) .aty. (y. with. grid) 

! OR 

pdf_at_y = EvalGridQucLnt(grid,pdf ( : , -6 : 6) ,y) 

with the latter being more eflicient if one needs to extract all flavours simultaneously. Note 
that here we have explicitly specifled the flavours, -6:6, that we wantJil 

5.1.2 Evolution representation 

For the purpose of carrying out convolutions, the human representation is not very advan- 
tageous because the splitting matrix in flavour space is quite complicated. Accordingly 
HOPPET uses a different representation of the flavour internally when carrying out con- 
volution of splitting matrices with PDFs. For most purposes the user need not be aware 
of this. The two exceptions are when a user plans to create derived splitting matrices 
(being careless about the flavour representation will lead to mistakes) or wishes to carry 
out repeated convolutions for a flxed nj value (appropriate manual changes of the flavour 
representation can speed things up). 

The splitting matrix can be simplifled considerably by switching to a different flavour 
representation, as can be seen in Table [H When carrying out a convolution, the only 
non-diagonal part is the block containing indices 0, 1. This representation is referred to 
as the evln representation. Whereas the human representation is rij-independent the evln 
depends on Uf through the S and g^g entries and the fact that flavours beyond n f are left 
in the human representation (since they are inactive for evolution with Uf flavours). 

To take a PDF in the human representation and make a copy in an evln representation, 
one uses the CopyHumanPdf ToEvln routine 

real(dp), pointer :: pdf _humaji( : , : ) , pdf _evln( : , : ) 

integer :: nf_lcl ! NB: nf would conflict with global variable 

[... setting up pdf_human, nf_lcl, etc. ...] 

call AllocPDF(grid,pdf _evln) ! or it might be an automatic array 

call CopyHumanPdf ToEvln(nf_lcl, pdf_humaji, pdf_evln) ! From human to evolution representation 

where one specifles the Uf value for the evln representation. One can go in the opposite 
direction with 

^If instead we had said pdf ( : , :) the resuh would have corresponded to a shce of flavours -6:7, i.e. 
including an interpolation of the representation labelling information, which would be meaningless. 
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call CopyEvlnPdf ToHuman(iif _lcl, pdf_evln, pdf_human) 

At any time one can check which is the representation of a given PDF using the GetPdf Rep 
function, 

integer nf_rep 

real (dp), pointer :: pdf ( : , : ) 

[ . . . set up pdf , . . . ] 
nf_rep = GetPdf Rep (pdf ) 

which returns the number of active flavours if the PDF is in the evln representation, or a 
negative integer if the PDF is in the human representation. 

5.2 Splitting function matrices 

Sphtting function matrices and their actions on PDFs are defined in module dglap_ob j ects 
(accessible as usual from module hoppet_vl). They have type split_mat. Below we shall 
discuss routines for creating specific predefined DGLAP splitting matrices, but for now we 
consider a general splitting matrix. 

The allocation of split_mat objects, 

type (split_mat) :: P 

integer : : nf _lcl 

call AllocSplitMat(grid, P, nf_lcl) 

is similar to that for grid_conv objects. The crucial difference is that one must supply 
a value for n/, so that when the splitting matrix acts on a PDF it knows which flavours 
are decoupled. From the point of view of subsequent initialisation a split_mat object just 
consists of a set of splitting functions. If need be, they can be initialised by hand, for 
example 

call InitGridConv(grid,P%qq , P_f unction_qq ) 

call InitGridConv(grid,P%qg , P_f unction_qg ) 

call InitGridConv(grid,P%gq , P_f unction_gq ) 

call InitGridConv(grid,P%gg , P_f unction_gg ) 

call InitGridConv(grid,P%NS_plus , P_f unction_NS_plus ) 
call InitGridConv(grid,P%NS_ininus, P_f unction_NS_minus) 
call InitGridConv(grid,P%NS_V , P_f unction_NS_V ) 

We can then write 

real(dp), pointer :: q(:,:), delta_q(:,:) 
[... allocations, etc. ...] 
delta_q = P .conv. q 
! OR 

delta_q = P * q 
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and delta_q will have the following components 

f5J:\ ^ ( P%qq P%qg \ /S\ 

5q+^ . = P%NS_plus ® g+g . (20) 
Sq^s,i = PyoNS_minus ® q^^ ^ 
6q^s = P°/.NS_V®g^g 

We have written the result in terms of components in the evln representation (and this 
is the representation used for the actual convolutions). When a convolution with a PDF 
in human representation is carried out, the program automatically copies the PDF to the 
evln representation, carries out the convolution and converts the result back to the human 
representation. The cost of changing a representation is O {Nx), whereas the convolution 
is O [N^), so in principle the former is negligible. In practice, especially when aiming for 
high speed at low A^^, the change of representation can imply a significant cost. In such 
cases, if multiple convolutions are to be carried out, it may be advantageous to manually 
change into the appropriate evln representation, carry out all the convolutions and then 
change back to the human manually representation at the end, see Sect. 15. 1[ 

As for grid_conv objects, a variety of routines have been implemented to help manip- 
ulate splitting matrices: 



type(split_mat) :: PA, PB, PC 
real (dp) : : factor 








call InitSplitMatCPA.PB [.factor] ) 


! PA 


= PB [*factor] 


(opt . alloc) 


call SetToZero(PA) 

call Multiply (PA, factor) 

call AddWithCoeff (PA, PB[, factor]) 


! PA 
! PA 
! PA 


= 

= PA * factor 
= PA + PB [*fact 


or] 


call SetToCoiivolutioii(PA,PB,PC) 
call SetToCommutator(PA,PB,PC) 


! PA 
! PA 


= PB*PC 

= PB*PC-PC*PB 


(opt . alloc) 
(opt . alloc) 


call Delete (split_mat) 


! PA 


s memory freed 





5.2.1 Derived splitting matrices 

As with grid_conv objects, HOPPET provides means to construct a split_mat object that 
corresponds to an arbitrary series of split_mat operations, as long as they all involve the 
same value of Uf. One proceeds in a very similar way as in Sect. 14.3.21 

real (dp), pointer :: probes( : , : , : ) 
type(split_mat) :: PA, PB, Pcomm 
integer : : i 

[. . .set nf_lcl, . . .] 
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call GetDerivedSplitMatProbes(grid,nf_lcl, probes) ! get the probes 

do i = 1 , size(probes,dim=3) ! carry out operations on each probe 

probesC : , : ,i) = PA* (PB*probes( : , : , i) ) - PB* (PA*probes( : , : , i) ) 
end do 

call AllocSplitMat(grid,Pcomm,nf _lcl) ! provide nf info in initialisation 

call SetDerivedConv(Pcomm, probes) ! Presult = [Pqg,Pgq] 

Note that we need to provide the number of active quark flavours to GetDerivedSplitMatProbes. 
As in Sect. 14.3.21 we first need to set up some 'probe' PDFs (note the extra dimension com- 
pared to earher, since we also have flavour information; the probe index always corresponds 
to the last dimension); then we act on those probes; finally we allocate the splitting matrix, 
and set its contents based on the probes, which are then automatically deallocated. 

5.3 The DGLAP convolution components 
5.3.1 QCD constants 

The splitting functions that we set up will depend on various QCD constants {uf, colour 
factors), so it is useful to here to summarise how they are dealt with within the program. 

The treatment of the QCD constants is not object oriented. There is a module (qcd) 
that provides access to commonly used constants in QCD: 

real(dp) :: ca, cf, tr, nf 
integer : : nf _int 

real (dp) :: betaO, betal, beta2 
[ ... ] 

Note that nf is in double precision 
To set the value of nj, call 

integer : : nf _lcl 
call qcd_SetNf (nf _lcl) 

where we have used the local variable nf _lcl to avoid conflicting with the nf variable 
provided by the qcd module. Whatever you do, do not simply modify the value of the nf 
variable by hand — when you call qcd_SetNf it adjusts a whole set of other constants (e.g. 
the P function coefficients) appropriately. 

There are situations in which it's of interest to vary the other colour factors of QCD, for 
example, if these colour factors are to be determined from a fit to deep-inelastic scattering 
experimental data. For that purpose, use 

real(dp) : : ca_lcl, cf_lcl, tr_lcl 
call qcd_SetGroup(ca_lcl, cf_lcl, tr_lcl) 

Again all other constants in the qcd module will be adjusted. A word of caution: the NNLO 
splitting functions actually depend on a colour structure that goes beyond the usual Ca, 
Cp and Tfj, namely dabcd"'^'^, which in the present version of hoppet is hard- wired to its 
default QCD value. 



— if you want the integer value of n/, use nf _int. 
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5.3.2 DGLAP splitting matrices 

The module dglap_objects includes a number of routines for providing access to the 
split_mat objects corresponding to DGLAP splitting functions 

type(split_mat) :: P_LD, P_NLD, P_NNLO 



! MSbar unpolarised case 

call InitSplitMatLD (grid, P_LO) 

call InitSplitMatNLO (grid, P_NLO) 

call IiiitSplitMatNNLO(grid, P_NNLO) 

! the MSbar polarised case. . . 

call InitSplitMatPolLO (grid, Pp_LD) 

call IiiitSplitMatPolNLO(grid, Pp_NLO) 

In each case the splitting function is set up for the nf and colour-factor values that are 
current in the qcd module, as set with the qcd_SetNf and qcd_SetGroup subroutine calls. 
If one subsequently resets the n/ or colour factor values, the split_mat objects continue 
to correspond to the Uf and colour factor values for which they were initially calculated. 
With the above subroutines for initialising DGLAP splitting functions, the normalisation 
is as given in eq. (jlj). 

In practice, because convolutions take a time O (A^), where A is the number of points 
in the grid, whereas additions and multiplications take a time 0{N), in a program it is 
more efficient to first sum the splitting matrices and then carry out the convolution, 

type(split_mat) : : P_sum 
real(dp), pointer :: q(:,:), dq(:,:) 



Note the use of brackets in the line setting dq: all scalar factors are first multiplied together 
[O (1)) so that we only have one multiplication of a PDF [O [N^)). Note also that we have 
chosen to include the (as2pi * dt) factor as multiplying the pdf, rather than the other 
option of multiplying P_sum, i.e. 

call Multiply(P_sum, (as2pi * dt)) 
dq = P_siim .conv. q 

The result would have been identical, but splitting matrices with positive interpolation 
order essentially amount to an O {7 x order x A) sized array, whereas the PDF is an 
O (13A) sized array and the for high positive orders that are sometimes used, it is cheaper 
to multiply the latter. 

The default for the NNLO splitting functions are the interpolated expressions, which 
are very fast to evaluate. Other possibilities, like the exact splitting functions or a previous 
set of approximated NNLO splitting functions which was used before the full calculation 



type(split_mat) :: Pp_LD, Pp_NLO 



! polarised 



call IiiitSplitMat(P_sum, P_LD) 
call AddWithCoeff (P_sum, P_NLD, as2pi) 
dq = (as2pi * dt) * (P_sum .conv. q) 
call Delete (P_sum) 



P_sum = P_LO 

P_sum = P_sum + as2pi * P_NLO 
Step dt in evolution 
Memory freed 
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was available are described in Appendix [Dl Note that the QCD colour factors introduced 
in Sect. 15.3.11 cannot be modified if the interpolated NNLO splitting functions are used, 
since these expressions use the default QCD values. 

5.3.3 Mass threshold matrices 

Still in the dglap_objects module, we have a type dedicated to crossing heavy quark mass 
thresholds. 

type(grid_d.ef ) : : grid 

type(mass_tliresliold_mat) : : MTM_NNLO 

call InitMTMNNLOCgrid, MTM_NNLQ) ! MTM_NNLO is coeff of (as/2pi)**2 

This is the coefficient of {as/lnY for the convolution matrix that accounts for crossing a 
heavy fiavour threshold in MS factorisation scheme, ai ixp = rrih, where ruh is the heavy- 
quark pole mass, as has been described in Sect. [21 Since the corresponding NLO term is 
zero, the number of fiavours in is immaterial at NNLO. 

The treatment of rij in the mass_threshold_mat is very specific because at NNLO, the 
only order in the MS factorisation scheme at which it's non-zero and currently known, 
it is independent of Uf. Its action does of course however depend on Uf. Since, as for 
splitjnat objects, we don't want to the action of of the mass_threshold_mat to depend 
on the availability of the current Uf information from the qcd module, instead we require 
that before using a mass_threshold_mat, you should explicitly indicate the number of 
fiavours (defined as including the new heavy fiavour). This is done using a call to the 
SetNfMTM(MTM_NNLO,nf_incl_heavy) subroutine. So for example to take a PDF in the 
effective theory with Uf = 3 active flavours pdf_nf3, and convert it to the corresponding 
PDF in the effective theory with Uf = 4 active flavours pdf_nf4 at m|, one uses code 
analogous to the following 

real(dp) :: pdf _nf 3( : , : ) , pdf_nf4(:,:) 
[ ... ] 

call SetNfMTM(MTM, 4) 

pdf_nf4 = pdf_nf3 + (as2pi)**2 * (MTM_NNLD. conv.pdf _nf3) 

The convolution only works if the pdf 's are in the human representation and an error is 
given if this is not the case. Any heavy flavour (like for example intrinsic charm) present 
in pdf _nf 3 would be left in unchanged. 

Note that the type mass_threshold_mat is not currently suitable for general changes 
of flavour-number. For example if you wish to carry out a change in the DIS scheme or at 
a scale fir ^ then you have to combine a series of different convolutions (essentially 
correcting with the lower number of flavours to the MS factorisation scheme at fip = rrih 
before changing the number of flavours and then correcting back to the original scheme 
and scale using the higher number of flavours) . 

As for the NNLO splitting functions, the mass threshold corrections come in exact and 
parametrised variants. By default it is the latter that is used (provided by Vogt [33]). The 
cost of initialising with the exact variants of the mass thresholds is much lower than for the 
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exact NNLO splitting functions (partly because there is no ri/ dependence, partly because 
it is only one flavour component of the mass-threshold function that is complex enough 
to warrant parametrisation). The variant can be chosen by the user before initialising the 
mass_threshold_mat by making the following subroutine call: 

integer : : threshold_variELnt 

call dglap_Set_nnlo_nf threshold(threshold_variaiit) 
with the following variants deflned (as integer parameters), again in the module dglap_choices: 
nnlo_nf threshold_exact 

nnlo_nf threshold_parain [default] 
5.3.4 Putting it together: dglap_holder 

The discussion so far in this subsection was intended to provide the reader with an overview 
of the different DGLAP components that have been implemented and of how they can be 
initialised individually. This is useful above all if the user needs to tune the program to 
some speciflc unusual application. 

In practice, one foresees that most users will need just a standard DGLAP evolution 
framework, and so will prefer not need to manage all these components individually. Ac- 
cordingly HOPPET provides a type, dglap_holder which holds all the components required 
for a given kind of evolution. To initialise all information for a flxed-flavour number evo- 
lution, one does as follows 

use hoppet_vl 

type(dglap_holder) : : dglap_h 

integer : : factscheme, nloop, nf_lcl 

nloop = 3 ! NNLO 

factscheme = f actscheme_MSbar ! or: f actscheme_DIS; f actscheme_PolMSbar 



! now do the initialisation 

call InitDglapHolderCgrid, dglap_h, factscheme, nloop) 

The constants f actscheme_* are deflned in module dglap_choices. The corrections to 
the splitting functions to get the DIS scheme are implemented by carrying out appropriate 
convolutions of the MS splitting and coefficient functions. Currently the DIS scheme is 
only implemented to NLOO. The polarised splitting functions are only currently known to 



Initialisation can also be carried out with a single call for a range of different numbers 
of flavours: 

®It's NNLO implementation would actually be fairly straightforward given the parametrisation provided 
in [34], and may be performed in future releases of hoppet. 



nf_lcl = 4 

call qcd_SetNf (nf _lcl) 
! call qcd_SetGroup( . . . ) 



! set the fixed number of flavours 

! if you want different colour factors 



NLO. 



26 



integer : : nflo, nfhi 
[. . .] 

nflo = 3; nfhi = 6 ! [calls to qcd_SetNf handled automatically] 
call InitDglapHolderCgrid, dglap_h, factscheme, nloop, nflo, nfhi) 

Mass thresholds are not currently correctly supported in the DIS scheme. 

For all the above calls, at NNLO the choice of exact of parametrised splitting func- 
tions and mass thresholds is determined by the calls to dglap_Set_nnlo_splitting and 
dglap_Set_nnlo_nf threshold, as described in Sects. 15.3.21 and 15.3.31 respectively. These 
calls must be made prior to the call to InitDglapHolder. 

Having initialised a dglap_holder one has access to various components: 

type dglap_holder 

type (split_mat) , pointer :: allP(l:nloop, nflo:nfhi) ! FFNS: nf lo=nf hi=nf _lcl 
type(split_mat) , pointer :: P_LO, P_NLO, P_NNLD 
type(mass_threshold_mat) : : MTM2 
logical : : MTM2_exists 

integer : : factscheme, nloop 

integer : : nf 

[ ... ] 
end type dglap holder 

Some just record information passed on initialisation, for example factscheme and nloop. 
Other parts are set up once and for all on initialisation, notably the allP matrix, which 
contains the 1-loop, 2-loop, etc. splitting matrices for the requested Uf range. 

Yet other parts of the dglap_holder type depend on nj. Before accessing these, one 
should first perform the following call: 

call SetNfDglapHolder(dglap_h, nf_lcl) 
This creates links: 

dglap_h%P_LO => dglap_h7.allP(l,nf_lcl) 
dglap_hy.P_NLO => dglap_h%allP(2,nf _lcl) 
dglap_hy.P_NNLO => dglap_h%allP(3,nf _lcl) 

for convenient named access to the various splitting matrices, and it also sets the global 
(qcd) Uf value (via a call to qcd_SetNf ) and where relevant updates the internal Uf value 
associated with MTM2 (via a call to SetNfMTM). 

As with other types that allocate memory for derived types, that memory can be freed 
via a call to the Delete subroutine, 

call Delete (dglap_h) 



6 DGLAP evolution 

So far we have described all the tools that are required to perform DGLAP convolutions of 
PDFs. In this section we describe how the different ingredients are put together to perform 
the actual evolution. 
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6.1 Running coupling 

Before carrying out any DGLAP evolutions, one first needs to set up a running_coupling 
object (defined in module qcd_coupling): 

type(ruiiiiiiig_coupling) : : coupling 

real(dp) :: alfas, Q, quark_masses(4 : 6) , muMat ch_mQuark 
integer : : nloop, fixnf 

[. . . set parameters . . .] 

call InitRunningCouplingCcoupling [, alfas] [, Q] [, nloop] [, fixnf]& 

& [, quark_masses] [, muMat ch_mQuark] ) 

As can be seen, many of the arguments are optional. Their default values are as follows: 

Q = 91.2_dp 

alfas = 0.118_dp ! Value of coupling at scale Q 
nloop = 2 

fixnf = [.not. present] 

! charm, bottom, top 

quark_masses(4:6) = (/ 1 .414213563_dp, 4.5_dp, 175.0_dp /) ! Heavy quark pole masses 
muMatch_mQuark = 1.0_dp 

The running coupling object is initialised so that at scale Q the coupling is equal to alfas. 
The running is carried out with the nloop /3-function. If the fixnf argument is present, 
then the number of flavours is kept fixed at that value. Otherwise flavour thresholds are 
implemented at scales 

muMatch_mQuark * quark_masses(4 : 6) 

where the quark masses are pole masses. The choice to use pole masses (and their partic- 
ular default values) is inspired by the benchmark runs [19] in which hoppet results were 
compared to those of Vogt's moment-space code QCD-Pegasus [5]. The default value of 
the QCD coupling is taken to be close to the present world average [55] . 

To access the coupling at some scale Q one uses the following function call: 

alfas = Value (coupling, Q [, fixnf]) 

This is the value of the coupling as obtained from the Runge-Kutta solution of the nloop 
version of eq. ([5]) (the numerical solution is actually carried out for l/a^), together with 
the appropriate mass thresholds. For typical values of as{Mz) the coupling is guaranteed 
to be reliably determined in the range 0.5 GeV < Q < 10^^ GeV. The values of the (3 
function coefficients used in the evolution correspond to those obtained with the values of 
the QCD colour factors that were in vigour at the moment of initialisation of the coupling. 

In the variable flavour-number case, the fixnf argument allows one to obtain the 
coupling for fixnf flavours even outside the natural range of scales for that number of 
flavours. This is only really intended to be used close to the natural range of scales, and 
can be slow if one goes far from that range (a warning message will be output). If one is 
interested in a coupling that (say) never has more than 5 active flavours, then rather than 
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using the f ixnf option in the Value subroutine, it is best to initiahse the couphng with a 
fictitious large value for the top mass. 

Often it is convenient to be able to enquire about the mass information embodied in a 
running_coupling. For example in the PDF evolution below, all information about the 
location of mass thresholds is obtained from the running_coupling type. 

The quark pole mass for flavour if Iv can be obtained with the call 

pole_mass = QuarkMass(coupling, iflv) 

The range of scales, Qlo < Q < Qhi for which iflv is the heaviest active flavour is obtained 
by the subroutine call 

call QRangeAtNf (coupling, iflv, Qlo, Qhi [, muM_mQ] ) 

The optional argument muMjiQ allows one to obtain the answer as if one had initialised the 
coupling with a different value of muMatch_mQuark than that actually used. One can also 
establish the number of active flavours, nf_active, at a given scale Q with the following 
function: 

nf_active = NfAtQ (coupling, Q [, Qlo, Qhi] [, muM_mQ] ) 

As well as returning the number of active flavours, it can also set Qlo and Qhi, which 
correspond to the range of scales in which the number of active flavours is unchanged. The 
optional muM_mQ argument has the same purpose as in the QRangeAtNf subroutine. The 
last of the enquiry functions allows one to obtain the range of number of flavours covered 
in this coupling, nflo < Uf < nfhi: 

call Nf Range(coupling, nflo, nfhi) 

Finally, as usual, once you no longer need a running_coupling object, you may free 
the memory associated with it using the Delete call: 

call Delete (coupling) 



6.2.1 Direct evolution 

We are now, at last, ready to evolve a multi-flavour PDF. This is done by breaking the 
evolution into steps, and for each one using a Runge-Kutta approximation for the solution 
of a first-order matrix differential equation. The steps are of uniform size in a variable u 
that satisfies the following approximate relation 



For a 1-loop running coupling one has u = (In In (5^/A)//3o, which is the variable that 
appears in analytical solutions to the 1-loop DGLAP equation. The step size in u, du, can 
be set with the following call 



6.2 DGLAP evolution 




(21) 
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real (dp) : : du = 0.1_dp ! or some smaller value 
call SetDef aultEvolutionDu(du) 

The error on the evolution from the finite step size should scale as (du)^. With the default 
value of du = 0.1, errors are typically somewhat smaller than 10~^ (see Sect. [9] for the 
detailed benchmarks). 

To actually carry out the evolution, one uses the following subroutine call: 



dglap_h 
coupling 
initial_pdf ( : , : ) 
Q_init, Q_end 
nloop 
untie nf 



type (dglap_holder) 
type (running_coupling) 
real (dp), pointer 
real (dp) 
integer 
integer 
[. . .] 

call EvolvePDF(dglap_h, initial_pdf , coupling, Q_init, Q_end & 

& [, muR_Q] [, nloop] [, untie_nf] [, du] ) 

which takes a PDF array pdf and uses the splitting matrices in dglap_h to evolve it from 
scale Q_init to scale Q_end. By default the renormalisation to factorisation scale ratio is 
muR_Q = 1.0 and the number of loops in the evolution is the same as was used for the 
running coupling (the nloop optional argument makes it possible to override this choice). 
Variable fiavour-number switching takes place at the pole masses (maybe one day this 
restriction will be lifted) as associated with the coupling. 

If the dglap_holder object dglap_h does not support the relevant number of loops 
or flavours, the program will give an error message and stop. With the untie_nf option 
you can request that the number of flavours in the evolution be 'untied' from that in the 
coupling in the regions where dglap_h does not support the number of flavours used in the 
coupling. Instead the closest number of flavours will be used0 

Mass thresholds (NNLO) are implemented as described in Sect. [2) 

Pdf„, = Pdf„^_i + I ' y ) (dglap_h%MTM2 .conv. pdf„^._J , (22a) 




pdf = pdf„^ - (dglap_h7oMTM2 .conv. pdf (22b) 



when crossing the threshold upwards and downwards, respectively. Note that the two 
operations are not perfect inverses of each other, because the number of flavours of the pdf 
used in the convolution differs in the two cases. The mismatch however is only of order 
(NNNNLO), i.e. well beyond currently known accuracies. 



^For example if dglapJi was initialised with Uf = 3... 5 while the coupling has rif = 3... 6, then 
variable flavour number evolution will be used up to ti/ = 5, but beyond the top mass the evolution will 
carry on with 5 flavours, while the coupling uses 6 flavours. There probably aren't too many good reasons 
for doing this (other than for examining how much it differs from a 'proper' procedure). 
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A general remark is that crossing a flavour threshold downwards will nearly always 
result in some (almost certainly physically spurious) intrinsic heavy-flavour being left over 
below threshold. 

6.2.2 Precomputed evolution and the evln_operator 

Each Runge-Kutta evolution step involves multiple evaluations of the derivative of the 
PDFs, and the evolution between two scales may be broken up into multiple Runge-Kutta 
steps. This amounts to a large number of convolutions. It can therefore be useful to create 
a single derived splitting matrix that is equivalent to the whole evolution between the two 
scales. 

A complication arises because evolutions often cross flavour thresholds, whereas a de- 
rived splitting matrix is only valid for flxed nj. Therefore a new type has to be created, 
evln_operator, which consists of a linked list of splitting and mass threshold matrices, 
breaking an evolution into a chain of interleaved flxed-flavour evolution steps and flavour 
changing steps. An evln_operator is created with a call that is almost identical to that 
used to evolve a PDF: 

type (evln_operator) :: evop 

real(dp), pointer :: pdf _iiiit( : , : ) , pdf _end( : , : ) 
[. . .] 

call InitEvliiDperator(dglap_h, evop, coupling, Q_init, Q_end & 

& [, muR_Q] [, nloop] [, untie_nf] [, du] ) 

It can then be applied to PDF in the same way that a normal split_mat would: 

pdf_end = evop * pdf_init ! assume both pdfs already allocated 

! OR (alternative form) 
pdf_end = evop .conv. pdf_init 

As usual the Delete subroutine can be used to clean up any memory associated with an 
evolution operator that is no longer needed. 

7 Tabulated PDFs 

The tools in the previous section are useful if one knows that one needs DGLAP evolution 
results at a small number of predetermined Q values. Often however one simply wishes to 
provide a PDF distribution at some initial scale and then subsequently be able to access 
it at arbitrary values of x and Q. For this purpose it is useful (and most efficient) to 
produce a table of the PDF as a function of Q^, which then allows for access to the PDF 
at arbitrary x, Q using an interpolation. All types and routines discussed in this section 
are in the pdf .tabulate module, or accessible also from hoppet_vl. 

7.1 Preparing a PDF table 

The type that contains a PDF table is pdf _t able. It first needs to be allocated, 
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typeCpdf _table) : : table 
[. . .] 



call AllocPdf Table (grid, table, Qmin, Qmax & 

& [, dlnlnQ ] [, lnlnQ_order ] [, f reeze_at_Qmin] ) 

where one specifies the range of Q values to be tabulated, from Qmin to Qmax, and optionally 
the interpolation step size dlnlnQ in the variable In In Q/ (0.1 GeV) (default dlnlnQ = 0.07, 
sufficient for 10~^ accuracy), the interpolation order lnlnQ_order, equal to 3 by default, 
and finally whether PDFs are to be frozen below Qmin, or instead set to zero (the default 
is freeze_at_Qmin= . false . , i.e. they are set to zero) 1^1 

By default a tabulation knows nothing about Uf thresholds, which means that in the 
neighbourhood of thresholds the tabulation would be attempting to interpolate a discontin- 
uous function (at NNLO). To attribute information about Uf thresholds to the tabulation, 
set them up first in a running coupling object and then transfer them: 

call AddNf Inf oToPdf TableCtable , coupling) 

When interpolating the table (see below), the set of Q values for the interpolation will be 
chosen to as to always have a common n/ value. Note that AddNf Inf oToPdf Table may 
only be called once for an allocated table: if you need to change the information about Uf 
thresholds. Delete the table, reallocate it and then reset the Uf information. This is not 
necessary if you just change the value of the coupling. 

Given an existing table, ref _table, a new table, new_table, can be allocated with 
identical properties (including any Uf information) as follows 

call AllocPdf Table(new_table , ref_table) 

All of the above routines can be used with 1-dimensional arrays of tables as well (in 
AllocPdf Table the reference table must always be a scalar). 

A table can be filled either from a routine that provides the PDFs as a function of x 
and Q, or alternatively by evolving a PDF at an initial scale. The former can be achieved 
with 

call FillPdfTable_LHAPDF(table, LHAsub) 

where LHAsub is any subroutine with the LHAPDF interface, i.e., as shown earlier in 
Sect. l5Xn 

To fill a table via an evolution from an initial scale, one uses 

type (pdf _table) : : table 

type (dglap_holder) :: dglap_h 

type (running_coupling) :: coupling 

real (dp), pointer :: initial_pdf ( : , : ) 

^°Notc that the spacing and interpolation in Q are treated independently of what's done in the PDF 
evolution. A reason for this is that in the PDF evolution one uses a variable related to the running 
coupling, which is similar to In In Q but whose precise details depend on the particular value of the coupling. 
Using this in the tabulation would have prevented one from having a tabulation disconnected from any 
coupling. Unfortunately du and dlnlnQ are not normalised equivalently — roughly speaking for Uf = A 
they correspond to the same spacing if dlnlnQ ~ 0.7du. 
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real (dp) : : QO 

integer : : nloop 

integer : : untie_nf 

[. . .] 

call EvolvePdfTableCtable, QO, initial_pdf , dglap_h, coupling & 
& [, muR_Q] [, nloop] [, untie_nf] ) 

which takes an the initial_pdf at scale QO, and evolves it across the whole range of Q 
values in the table, using the EvolvePDF routine. The arguments have the same meaning 
as corresponding ones in EvolvePDF, explained in Sect. l6.2H The du value that's used is the 
default one for EvolvePDF which, we recall, may be set using SetDef aultEvolutionDu(du) . 
If the Q spacing in the tabulation is such that steps in du would be too large, then the 
steps are automatically resized to the tabulation spacing. 

One may also use the precomputed evolution facilities of Sect. I6.2.2[ by calling the 
routine 

call PreEvolvePdf Table(table , QO, dglap_h, coupling, & 

& [, muR_Q] [, nloop] [, untie_nf] ) 

prepares evln_operators for all successive Q intervals in the table. An accelerated evolu- 
tion that uses these operators instead of explicit Runge-Kutta steps may then be obtained 
by calling 

call EvolvePdfTableCtable, initial_pdf) 

The EvolvePdf Table routine may be called as many times as one likes, for different ini- 
tial PDFs for example; however, if one wishes to change the parameters of the evolution 
(coupling, perturbative order, etc.) in the precomputed option, one must first Delete the 
table and then prepare again it. 



7.2 Accessing a table 



The main way to access a table is as follows 

real(dp) :: pdf(-6:6), y, x, Q 
[. . .] 

call EvalPdfTable_yQ(table, y, Q, pdf) 
! or using x 

call EvalPdfTable_xQ(table, x, Q, pdf) 

There may be situations where it is useful to access the internals of a pdf _table, for 
example because one would like to carry out a convolution systematically on the whole 
contents of the table. Among the main elements are 



type pdf_table 
integer 

real (dp), pointer 
real (dp), pointer 
integer, pointer 
real (dp), pointer 
[. . .] 

end type pdf_table 



nQ 

tab(: , : , :) 
Q_vals(:) 
nf _int( : ) 
as2pi( : ) 



arrays run from 0:nQ 
the actual tabulation 
the Q values 
nf values at each Q 
alphas (Q)/2pi at each Q 
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where the third dimension of tab spans the range of tabulated Q values, and the Uf and 
coupling information are only allocated and set if one has called AddNflnfoToPdf Table 
for the table. 

An example of usage of the low-level information contained in the table is the following, 
which initialises table_deriv_LO with the LO derivative of table: 

do iQ = 0, table7.nq 

table_deriv_L07.tab( : , : ,iQ) = table7.as2pi(iC)) * & 

& ( dglap_h°/.allP(l,table%nf _int(iQ)) * table%tab( : , : , iQ) ) 

end do 

where we assume table_deriv_LO to have been allocated with an appropriate structure at 
some point, e.g. via 

call AllocPdfTable(table_deriv_LD, table) 

The above mechanism has found use in the a-posteriori PDF library [T71 [T8] and in work 
matching event shapes with fixed-order calculations [151 [T3]. Oiis could also imagine using 
it to obtain tables of (flavour-separated) structure functions, if one were to convolute with 
coefficient functions rather than splitting functions. 

As with all other objects, a pdf _table object can be deleted using 

call Delete (table) 

Notice that a table is a local variable in each procedure and so effectively a different variable 
each separate procedure. If ones needs to use a PDF table across different procedures, it 
has to be defined within a module. In Appendix \^ there is a detailed example of different 
ways of accessing a table. 

8 Streamlined interface 

Now we present the streamlined interface to hoppet, intended to allow easy access to the 
essential evolution functionality from languages other than F95. It hides all the object- 
oriented nature of the program, and provides access to one pdf _table, based on a single 
grid definition. The description will be given as if one is calling from F77. An include file 
src/hoppet_vl .h is provided for calling from C++ — the interface is essentially identical 
to the Fortran one, with the caveat that names are case sensitive (the cases are as given 
below), and that PDF components referred to below as pdf (-6:6) become pdf [0. .12]. 
A summary of the most relevant procedures of this interface and their description can be 
found in the reference guide. Appendix [Bl 

8.1 Initialisation 

The simplest way of initialising the streamlined interface is by calling 

call hoppetStartCdy.nloop) 
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which will set up a compound (four different spacings) grid with spacing dy at small x, 
extending to y = 12, and numerical order = — 5. The Q range for the tabulation will be 
1 GeV < Q < 28 TeV and a reasonable choice will be made for the dlnlnQ spacing (related 
to dy). It will initialise splitting functions up to nloop loops (though one can still carry out 
evolutions with fewer loops). If you need more control over the initialisation, you should 
use 

call hoppetStartExtendedCymax , dy , Qmin , Qmax , dlnlnQ , nloop , order , fact scheme) 

which will again set up compound grid, but give control over the numerical order and the 
y and Q ranges and spacings (as before dy is the spacing at small x). It also allows one to 
choose the type of evolution according to f actscheme. 

8.2 Usage 

To carry out an evolution, one should first decide whether one wants a fixed-flavour number 
scheme or a variable flavour number scheme (the default). Either can be set with its 
parameters as follows: 

call hoppetSetFFN(f ixed_nf ) 

call hoppetSetVFN(mc , mb, mt) ! Heavy quark pole masses 

where for the VFN one specifies the pole-masses for the quarks. An evolution is carried 
out with the following routine 

call hoppetEvolve(asQ , QOalphas, nloop, muR_Q, LHAsub, QOpdf) 

where one specifies the coupling asQ at a scale QOalphas, the number of loops for the 
evolution, nloop, the ratio of renormalisation to factorisation scales muR_Q0 the name of 
a subroutine LHAsub with interface 

subroutine LHAsub (x,Q,pdf) 
implicit none 

double precision x,Q ,pdf (-6 : 6) 

[...] ! sets pdf to be momentum densities, e.g. pdf(O) = xg(x) 
end subroutine 

to return the initial condition for the evolution and the scale QOpdf at which one starts 
the PDF evolution. Note that the LHAsub subroutine will only be called with Q = QOpdf. 
To access the coupling one uses 

alphas = hoppetAlphaS(Q) 

^^Note that in the streamHned interface, with muR_Q 7^ 1, the running coupling flavour thresholds are 
still placed at the quark masses; the evolution needs the coupling for a given number of flavours outside the 
standard range for that number of flavours (precisely because muR_Q 7^ 1) and this is done automatically 
in the evolution. In contrast in the benchmark studies |19| . the flavour thresholds for the coupling were 
placed at muR_Q x mg. This is a perfectly valid alternative, but can complicate the specification of the 
as value — for example with muR_Q = 0.5 the matching for the top threshold would be carried out at 
fi = 0.5mt ~ 85 GeV, and if one specified the coupling at scale Mz, it wouldn't be clearer whether this 
was a 5-flavour value or a 6-flavour value. With the procedure adopted in the streamlined interface the 
issue does not arise. (While in F95 the user has the freedom to do as they prefer). 
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while the PDF at a given value of x and Q is obtained with 

call hoppetEval(x,Q,f ) 

which sets f (-6:6) (recall that it is xg{x), etc., that are returned, since this is what is 
used through hoppet). 

It is also possible to prepare an evolution in cached form. This is useful if one needs to 
evolve many different PDF sets with the same evolution properties (coupling, initial scale, 
etc.), as is the usual situation in global analyses of PDFs, because though the preparation 
may take a bit longer than a normal evolution (2-10 times depending on the order), once 
it is done, cached evolutions run 3-4 faster than a normal evolution. The preparation of 
the cache is carried out with 

call hoppetPreEvolveCasQ , QOalphas, nloop, muR_Q, QOpdf) 

and then the cached evolution is carried out with 

call lioppetCacliedEvolve(LHAsub) 

The results may be very slightly different from those in a normal evolution (some infor- 
mation is lost when caching), and the user may wish to check on a case-by-case basis that 
such differences don't matter in practice. 

The tabulation can also be filled with the contents of an external PDF package, by 
calling 

call lioppetAssigii(LHAsub) 

where LHAsub is the name of any subroutine with the interface given above, which will now 
be called with a range of Q values corresponding to the internal tabulation scales. This 
essentially just transfers an external tabulation into hoppet's internal representation. 

Finally given an evolved or assigned PDF, one can obtain information about convolu- 
tions of the splitting functions with the PDFs: 

call hoppetEvalSplitCx , Q , iloop.nf ,f ) 

sets f (-6 : 6) equal to the value at x, Q of convolution of the iloop splitting function matrix 
(with nf flavours) with the currently tabulated PDF. If nf < the number of flavours used 
is the one appropriate at the specified Q scale (as long as the information is available, i.e. 
one of hoppetEvolve or hoppetCachedEvolve has been called). The first call with a given 
nf for a specified iloop will be slow (~ the time for a cached evolution), but subsequent 
calls with the same values will be fast. 

The routines described here are to be found in src/streamlined_interf ace . f 90 and 
may provide inspiration for the user wishing to write their own F95 code for hoppet. 

9 Benchmarks 

Key questions in assessing the usefulness of a PDF evolution code include that of its 
correctness, its accuracy and its speed, hoppet's correctness has been established with a 
reasonable degree of confidence in the benchmark tests [19] where it was compared with 
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the Mellin space based evolution code QCD-Pegasus [3]. The program used to carry out 
those tests is available as benchmarks/benchmarks . f 90. The user should carefully read 
the detailed comments at the beginning for usage instructions. 

The results used in [19] were obtained with very finely spaced grids, in order to guar- 
antee small numerical errors (< 10"''). Such accuracies are useful when comparing and 
testing two independent codes, because differences or bugs in the implementation of the 
physics (especially the higher-order parts) may only manifest themselves as small changes 
in the results. 

In contrast, for use in most physical applications, an accuracy in the range 10"^ to 
10"^ is generally y more than adequate, since it is rare for other sources of numerical 
uncertainty (e.g. Monte Carlo integration errors in NLO codes, or experimental errors) to 
be comparably small. The critical issue in such cases is more likely to be the speed of the 
code, for example in PDF fitting applications. 

HOPPEt's accuracy and speed both depend on the choice of grid (in y) and the evolution 
and/or tabulation steps in Q. We shall start with the question of the accuracy. 



9.1 Accuracy 

To measure the accuracy, we use the same initial condition and evolution parameters as 
in [191: 
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xu^{x) = 5.107200x°-^(l - xf , (23a) 

xd^{x) = 3.064320x°-^(l - x)^ , (23b) 

xd{x) = 0.1939875x-°-^(l - xf , (23c) 

xu{x) = xd{x){l — x) , (23d) 

xs{x) = xs{x) = 0.2{xd{x) + xu{x)) , (23e) 

xg{x) = 1.7x-°-\l - xf , (23f) 

where Uy = u — u, dy = d — d, and all other flavours are zero. The initial scale i 
Qo = (v^ — e) GeV, as{Qo) = 0.35 and the charm, bottom and top pole masses are kept 
at the default values, V2, 4.5 and 175 GeV respectively (as used also in the streamlined 
interface). Observe that the initial conditions and coupling are actually both given for 
three active flavours (i.e. infinitesimally below the charm mass). The evolution is carried 
out to NNLO accuracy in a variable flavour number scheme, including the mass thresholds 
in the coupling and PDF. 

All tests here are carried out based on tabulations of the PDF evolution. Sect. [71 We 
first run HOPPET with a very fine grid spacing to provide a reference result. Then we 
run the evolution for a coarser grid — the accuracy of the coarser grid is determined 
by comparing its results with those from the reference grid. We determine the relative 
accuracy for each flavour at 5000 points in the x, Q plane, as shown in Fig. [2l The points 



""^^With e an infinitesimal number. Note that this is unrelated to the evolution accuracy e, introduced 
later in this section. 
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X 

Figure 2: The set of points in x, Q used to determine the accuracy of the evolution. The 
areas shaded in grey are regions where one of the flavours is in the neighbourhood of a 
sign-change (bottom-left: c, c, right: u) and so is ignored in the accuracy determination. 
The accuracies shown here correspond to a y grid with a base spacing of dy = 0.2 and 
other parameters as described in the text for Fig. [3l The colour coding indicates the error 
in the least-well determined (non-excluded) flavour channel (b . . .b) at each point. 

are uniformly spaced in C = lnl/x + 9(l — x) so as to obtain fine coverage at small and large 
X. The Q values are chosen more closely spaced at low Q where the evolution is fastest 
and they are taken slightly correlated with ( so as to cover a nearly continuous range of 

qE 

One difficulty that arises when examining relative accuracies is that some flavours 
change sign as one varies x oi Q. Close to the zero the relative accuracy diverges because 
of the small value of the denominator. Therefore in global accuracy estimates, we eliminate 
flavours in the region where they change sign (within A( = 0.4 and at the neighbouring Q 
value). Specifically, for our initial conditions, this corresponds to c, c for the two lowest Q 
values below x ~ 10~^ and u for x > 0.90 The exact regions are shaded in grey in Fig. [21 

^•^This procedure differs from that in [19] where fewer (500) points were used and one compared not 
individual flavours, but combinations intended to be more directly revealing of any deficiencies in the 
evolution. This reflects the difference in needs between obtaining a global measure of the accuracy and 
providing benchmarks intended in part to facilitate the debugging of independent codes. 

^^Thc change in sign of the charm distribution is not worrying physically since it is close to threshold 
where it will be compensated for by finite mass effects in the coefficient functions; for u the sign change 
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Our tabulation covers the range 10-5 < X < 1, < g < 10^ GeV. The Errid in 
y = In 1/x will consist of 4 nested subgrids: one covering the whole y range with spacing 
dy and others with spacings dy/3, dy/9, dy/27 extending to y = 2,0.5,0.2 respectively. 
Except where stated we shall use order = —6. In Q the default interpolation order will 
be 4. The reference grids use dy = 0.025 and dlnlnQ = 0.005. 

Fig. [3] shows the relative accuracy e as a function of x for two grid-spacing choices (left 
dy = 0.2, right dy = 0.05, dlnlnQ = dy/4 in both cases). Each solid line corresponds to 
one Q value and shows the error in the least- well-determined flavour channel at each x, 
excluding flavour channels close to a sign-change. The relative accuracy e is poorest as one 
approaches x = 1, where the PDFs all go to zero very rapidly and so have divergent loga- 
rithmic derivatives in x, dlnq/dlnx, adversely affecting the accuracy of the convolutions. 
This region is always the most difficult in x-space methods, however the use of multiple 
subgrids in x allows to one to obtain acceptable results for x < 0.9 which is likely to be 
the largest value of any phenomenological relevance. 

At X ~ 0.1, 0.6 and 0.8 one notices step-like structures — these are the points where one 
switches between subgrids, with a significant degradation in accuracy at x values below the 
transition. These structures are also visible in the colour-coded accuracy representation 
in Fig. [2], which corresponds to dy = 0.2 and allows one to visualise more clearly the Q 
dependence of the accuracy. The effect of the grid spacing is clearly visible as one goes 
from the left to the right-hand plots of Fig. [31 with the reduction in the spacings by a 
factor of 4 leading to an improvement in accuracy by a factor ~ 100. 

For completeness we also show the parts of the charm channel that have been excluded 
because of the proximity to a sign change (dashed lines, lower-left shaded region of Figj2]). 
One observes in particular a spike near x ~ 7 x 10"^ where the charm distribution has its 
zero. Including this in a estimate of the global accuracy would be senseless since it actually 
corresponds to a divergence and the peak- value for the spike is arbitrary, depending on the 
precise choice of points used to estimate the accuracy. The question of the exact region to 
exclude is somewhat arbitrary, but the choice made above seems not unreasonable in the 
light of Fig. H 

Fig. [3] is useful in order to obtain a detailed picture of the accuracy of the evolution 
with a given set of parameters. To quote a single, global, number for the relative accuracy 
e we make the conservative choice of taking the largest value of e that occurs in a chosen 
X range. We will examine a restricted range, x < 0.7, studying just the g,u, d, s flavours, 
and also a wider range, x < 0.9 with all flavours. 

Fig. m shows the effect of varying the base dy and dlnlnQ separately, while the other is 
fixed at the reference value. The 'guds' flavours in the x < 0.7 range are generally better 
determined, for a given set of grid parameters, than the full set of flavours up to x < 0.9. 
This is as one would expect since the large-x region is usually the hardest and the 'guds' 
flavours are generally somewhat smoother than the the others. Using only three y subgrids 

is more surprising, though it may be related to non-trivial interactions between the evolutions of the u 
and u components at NNLO (note that in the region of the sign change they differ by many orders of 
magnitude). 
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Figure 3: The relative accuracy e of the least well determined flavour channel at each a:, Q 
point, shown as a function of x for many Q values. The results for the part of the charm 
distribution excluded from the analysis (near sign change) are shown separately. 




4 subgrids (1 :3:9:27) 
Ssubgrids (1:3:15) 
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all, x<0.9 [o = 3] 
all, x<0.9 [0 = 4] 
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Figure 4: Left: the (globally) worst relative accuracy e as a function of the base y-grid 
resolution parameter, dy — shown for two 2/-grid configurations and two x-ranges and 
flavour-sets. Right: the relative accuracy as a function of the resolution in InlnQ of the 
tabulation, dlnlnQ, for different x/fiavour ranges and for different lnlnQ_order values (o). 



40 





lf95 


ifort 


g95 


ts 




0.9 


0.66 


2.8 


ta 


[ms] 


0.16 


0.12 


0.13 


ti 


[ms] 


37 


38 


330 


tp 


[ms] 


51 


44 


310 


tc 


[ms] 


8.8 


9.8 


110 


txQ 


[/is] 


2.7 


3.1 


25 



Table 2: Contributions to the run time in eqs. ( l24l) for dy = 0.2 and dlnlnQ = 0.05 and 
standard values for the other parameters (on a 3.4GHz Pentium IV (D) with 2 MB cache). 

worsens the situation when including the largest x values, and reducing the order in the Q 
interpolation also adversely affects the accuracy (also for x < 0.7, not shown). 

From Fig. HI we deduce that inaccuracies from the Q and y parts of the grid are similar 
when dlnlnQ = dy/4. This is the combination that we shall use as standard. 

9.2 Timing 

The time spent in hoppet for a given analysis can expressed as follows, according to 
whether or not one carries out pre-evolution: 

^no pre-ev tg -\- TT-ota -|- Tli(ti -\- TlxQ txo) , (24a) 
^with pre-ev = tg + na{ta + tp) + nj(tc + UxQ t^o) , (24b) 

where tg is the time for setting up the splitting functions, rta is the number of different 
running couplings that one has, t^ is the time for initialising the coupling, rij is the number 
of PDF initial conditions that one wishes to consider, ti is the time to carry out the 
tabulation for a single initial condition, is the number of points in x, Q at which one 
evaluates the full set of flavours once per PDF initial condition; in the case with pre- 
prepared cached evolution, tp is the time for a preparing a cached evolution and t^. is the 
time for performing the cached evolution. Finally t^g is the time it takes to evaluate the 
PDFs at a given value of (x, Q^) once the tabulation has been performed. 

The various contributions to the run-time are shown in table [2] for dy = 0.2 and 
dlnlnQ = 0.05 (giving an accuracy ~ 10""^), for various compilers. In a typical analysis 
where run-times matter, such as a PDF fit, it is to be expected that the time will be 
dominated by tc (or ti). However, in the typical case of global PDF fits, for which the 
number of x,Q points is rather large (> 3000), it will be n^Qt^Q that takes the most 
timj^. We note (with regret!) the considerably speed advantage (almost an order of 
magnitude) that is to be had with commercial compilers. 

We study t^ and ti in more detail in Fig. where we relate them to the accuracy 
obtained from the evolution. As one would expect, studying just the 'guds' flavours for 

With these numbers, it is easy to check that a global fit with n^Q ~ 3000 and rii ~ 10^ could be 
completed in less than half and hour. 
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Figure 5: relative accuracy obtained as a function of the time taken to perform the evolution 
(lf95, 3.4GHz Pentium IV (D) with 2 MB cache). 



X < 0.7 one obtains better accuracy for a given speed than with all flavours for x < 0.9. 
Overall one can obtain 10~^ accuracy with tc — 10"^ s and 10~^ accuracy with — 10"^ s. 
In general ti is about 4—5 larger than tc, highlighting the advantage of the cached evolution. 

We note that the time t^Q for evaluating each point is essentially independent of dy 
and dlnlnQ. If a computation is dominated by t^Q., then it can be made somewhat faster 
by lowering the interpolation orders, at the expense of needing a finer grid (and so longer 
evolution times). 

The timings shown here are roughly similar, for accuracies ~ 10~^, to those obtained 
with the A^-space code Pegasus [3] when the number of x, Q points to be evaluated is 
(9(10^). For much smaller numbers of points Pegasus becomes superior (because of the 
significantly smaller ratio tc/ixq), while for much larger numbers of points hoppet becomes 
better. Other NNLO evolution codes published in recent years [5], El [T] are generally less 
competitive either in terms of accuracy or speed. 

To close this section, we summarise in table [3] the different parameters that are relevant 
in determining the accuracy of the evolution and tabulation, together with comments about 
the components of the timing affected by each parameter. 



10 Conclusions 

HOPPET is an x-space evolution code that is novel both in terms of the accuracy and speed 
that it provides compared to other x-space codes, and in terms of its interface, designed 
to provide a straightforward and physical way of manipulating PDFs beyond the built-in 
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Parameter 


Default 


Timing impact 


Notes 


base dy 
order 

Def aultConvolutionEps 


[-6] 


tsi til tpi tc 
ts: tp: txQiti: tc) 
ts 


Default subgrids in ratio 1:3:9:27 
final acc. limited by ~ twice this. 


du 

dlnlnQ 
lnlnQ_order 


0.1 

[dy/4] 

4 


ti ) tp 
til tpi tc 
txQ 


immaterial if > 1.4 dlnlnQ 


Def aultCouplingDt 


0.2 


ta 


default sufficient for e ~ 10~^ 



Table 3: Parameters involved in the accuracy of a tabulated evolution. Default values 
shown in square brackets are to be specified by hand in the F95 interface, but are auto- 
matically set in the simpler of the initialisation calls with the streamlined interface. 



task of DGLAP evolution. 

Features that might be envisaged for future releases include DIS coefficient functions, 
full support for the DIS factorisation scheme, and the addition of time-like evolution, 
relevant for phenomenological fits to fragmentation functions as in [36]. In principle, the 
information presented here is sufficient to allow a user to implement the coefficient functions 
themselves, while the DIS scheme and timelike evolution would require somewhat more 
knowledge of the internals of the program. 

More ambitious possible extensions cover a wide range of physics. Just within QCD, a 
general physical feature absent from mainstream PDF evolution codes is that of evolution 
that includes matching with various types of resummed calculations. Although studies 
in this direction have already been performed, both for small x resummations, as in [M] 
and for large x resummation, as in [37], no general public code exists which performs a 
matching between resummed and fixed (NLO, NNLO) order splitting functions, either in 
the time-like case or in the space-like case. 

Another interesting extension of HOPPET would be to implement a more general mass 
treatment of heavy quarks. A proper treatment of heavy quark mass effects is required 
to obtain a good description of heavy flavour structure function as measured in HERA. 
Although there are by now several studies of the effect of heavy quark masses in global 
analysis of PDFs [38l 139]. which show sizable effects on predictions for LHC observables, 
there does not exist right now a public code were this General Mass heavy quark schemes 
are implemented. 

Also of interest are non-QCD effects. Evolutions with QED radiation have been pre- 
sented in [5l 130], however so far no public code exists for evolution including both QCD 
and electroweak (EW) effects [ITl 1^ . This is of particular relevance at LHC energies, 
since flavour is associated with an SU(2) charge and so soft divergences (above M^) do 
not cancel in the PDF evolution between real and virtual contributions, leading one to 
expect non-negligible effects in the flavour structure of the PDFs at high scales. The full 
QCD+EW evolution is a rather task complicated because of the need to include the EW 
flavour couplings, including the CKM matrix, and polarisation. For this kind of problem 
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a code such as hoppet provides a good starting point, since it has a clean separation 
of the numerical and flavour aspects of evolution, and verified unpolarised and polarised 
evolution. 
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A Example programs 
A.l General interface 

The program below generates a subset of table 15 of the NNLO benchmark evolution 
in the second reference of jT9]. It is to be found (in a slightly more commented form) 
in example_f 90/tabulation_exainple . f 90. Compilation instructions are to be found 
in the README file in the main directory of the release. A program that has the same 
functionality but in F77, using the streamlined interface of Sect. [HI is to be found as 
example_f 77/tabulation_example . f . 

progremi tabulatioii_exajiiple 
use hoppet_vl 
implicit none 

real(dp) :: dy, ymax, quark_masses(4 : 6) 

integer : : order, nloop, ix 

type (grid_def ) :: grid, gdarray(4) ! holds information about the grid 
type (dglap_holder) : : dh ! holds the splitting functions 

type (pdf _table) :: table ! holds the PDF tabulation 

type(running_coupling) : : coupling 

real (dp), pointer :: pdfO(:,:) ! holds the initial pdf 

real(dp) :: QO, Q, pdf _at_xq(-6 : 6) 

real(dp), parameter :: heralhc_xvals(9) = & 

& (/le-5_dp , le-4_dp , le-3_dp , le-2_dp , . l_dp , . 3_dp , . 5_dp , . 7_dp , . 9_dp/) 

! set up parameters for grid 
order = -6 
ymax = 12.0_dp 
dy = 0. l_dp 

! set up the grid itself — we use 4 nested subgrids 

call InitGridDef (gdarray(4) ,dy/27.0_dp, 0.2_dp, order=order) 

call InitGridDef (gdarrayO) ,dy/9.0_dp, 0.5_dp, order=order) 
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call InitGridDef (gdarray(2) ,dy/3.0_dp, 2.0_dp, order=order) 
call InitGridDef (gdarray(l) ,dy, ymax , order=order) 

call InitGridDef (gr id , gdarray (1:4), locked= . true . ) 

! initialise the splitting-function holder 
nloop = 3 

call InitDglapHolder (grid , dh , f actscheme=f actscheme_MSbar , & 
& nloop=nloop,nf lo=3 ,nf hi=6) 

! initialise a PDF from the function below (must be contained, 
! in a "used" module, or with an explicitly defined interface) 
call AllocPDF(grid, pdfO) 

pdfO = unpolarized_dummy_pdf (xValues(grid) ) 
QO = sqrt(2.0_dp) ! the initial scale 

! allocate and initialise the running coupling with a given 
! set of quark masses (NB: charm mass just above QO) . 
quark_masses(4:6) = (/I .414213563_dp, 4.5_dp, 175.0_dp/) 
call InitRunningCoupling(coupling , alf as=0 . 35_dp , Q=QO , nloop=nloop , & 
& quark_masses = quark_masses) 

! create the tables that will contain our copy of the user's pdf 
! as well as the convolutions with the pdf. 

call AllocPdfTable(grid, table, Qmin=1.0_dp, Qmax=10000.0_dp, & 

& dlnlnQ = dy/4.0_dp, f reeze_at_Qmin= . true . ) 
! add information about the nf transitions to the table (improves 
! interpolation quality) 
call AddNf Inf oToPdf Table(table , coupling) 

! create the tabulation based on the evolution of pdfO from scale QO 

call EvolvePdfTable(table, QO, pdfO, dh, coupling, nloop=nloop) 

! alternatively "pre-evolve" so that subsequent evolutions are faster 

! call PreEvolvePdf Table(table, QO, dh, coupling) 

!call EvolvePdf Table (table, pdf 0) 

! get the value of the tabulation at some point 
Q = 100.0_dp 

write(6, ' (a,f8.3,a) ') " Evaluating PDFs at Q = ",Q," GeV" 

write(6, ' (a5,2al2,al4,al0,al2) ' ) "x",& 

& "u-ubar" , "d-dbar" , "2 (ubr+dbr) " , "c+cbar" , "gluon" 
do ix = 1, size (heralhc_xvals) 

call EvalPdf Table_xQ(table,heralhc_xvals(ix) , Q ,pdf _at_xQ) 
write(6, ' (es7. I,5esl2.4) ') heralhc_xvals(ix) , & 

& pdf _at_xQ(2)-pdf _at_xQ(-2) , pdf _at_xQ(l)-pdf _at_xQ(-l) , & 

& 2*(pdf_at_xQ(-l)+pdf_at_xQ(-2)) , (pdf _at_xQ(-4)+pdf _at_xQ(4) ) , & 

& pdf_at_xQ(0) 

end do 

! some clecining up (not strictly speaking needed, but illustrative) 
call Delete (table) ; call Delete(pdf 0) ; call Delete(dh) 
call Delete (coupling) ; call Delete (grid) 
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contains 



! ! The dummy PDF suggested by Vogt as the initial condition for the 
! ! unpolarized evolution (as used in hep-ph/0511119) . 
function unpolarized_dummy_pdf (xvals) result (pdf) 
real (dp), intent (in) :: xvals(:) 

real(dp) :: pdf (size (xvals) , -6 : 7) ! note upper bound! 

real(dp) :: uv (size (xvals) ) , dv(size(xvals) ) 
real(dp) :: ubar(size(xvals) ) , dbar (size (xvals) ) 



real (dp), parameter :: N_g = 1.7_dp, N_ls = 0.387975_dp 
real (dp), parameter :: N_uv=5 . 107200_dp, N_dv = 3.064320_dp 
real (dp), parameter :: N_db = half*N_ls 



pdf = zero 



! — remember that these are all xvals*q(xvals) 
uv = N_uv * xvals**0 . 8_dp * (l-xvals)**3 
dv = N_dv * xvals**0.8_dp * (l-xvals)**4 
dbar = N_db * xvals** (-0 . l_dp) * (l-xvals)**6 
ubar = dbar * (1-xvals) 



! labels iflv_g, etc., come from the hoppet_vl module 



pdf ( 


, iflv. 


-g) 


= N_g * xvals** (- 


-0.1_dp) 


pdf ( 


,-iflv. 


.s) 


= 0.2_dp*(dbar + 


ubar) 


pdf ( 


, iflv. 


.s) 


= pdf ( : , -if lv_s) 




pdf ( 


, iflv. 


.u) 


= uv + ubar 




pdf ( 


,-iflv. 


.u) 


= ubar 




pdf ( 


, iflv. 


.d) 


= dv + dbar 




pdf ( 


,-iflv. 


.d) 


= dbar 





end function unpolarized_dummy_pdf 



end progrsun tabulation_example 

The expected output from the program is: 











Evaluatin 


g PDFs at Q 




100.000 GeV 










X 




u-ubar 




d-dbar 


2(ubr+dbr) 




c+cbar 




gluon 


1 


OE- 


05 


3 


1907E-03 


1 


.9532E-03 


3 


4732E+01 


1 


5875E+01 


2 


2012E+02 


1 


OE- 


04 


1 


4023E-02 


8 


.2749E-03 


1 


5617E+01 


6 


7244E+00 


8 


8804E+01 


1 


OE- 


03 


6 


0019E-02 


3 


.4519E-02 


6 


4173E+00 


2 


4494E+00 


3 


0404E+01 


1 


OE- 


02 


2 


3244E-01 


1 


.3000E-01 


2 


2778E+00 


6 


6746E-01 


7 


7912E+00 


1 


OE- 


01 


5 


4993E-01 


2 


.7035E-01 


3 


8526E-01 


6 


4466E-02 


8 


5266E-01 


3 


OE- 


01 


3 


4622E-01 


1 


.2833E-01 


3 


4600E-02 


4 


0134E-03 


7 


8898E-02 


5 


OE- 


01 


1 


1868E-01 


3 


.0811E-02 


2 


3198E-03 


2 


3752E-04 


7 


6398E-03 


7 


OE- 


01 


1 


9486E-02 


2 


.9901E-03 


5 


2352E-05 


5 


6038E-06 


3 


7080E-04 


9 


OE- 


01 


3 


3522E-04 


1 


.6933E-05 


2 


5735E-08 


4 


3368E-09 


1 


1721E-06 



The file example_f 90/tabulation_example .def ault_output contains a copy of these re- 
suhs, so as to allow easy comparison. The numbers correspond to evolution with variable 
flavour number, hf = /^r, and the parametrised versions of the NNLO splitting functions 
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and mass threshold terms. The reader may verify that they correspond to those given in 
the top panel of table 15 of the second reference of [19] {^f = f^R,)- 



A. 2 Streamlined interface 

The program below generates exactly the same output as the previous example program, 
but this time using the streamlined interface introduced in Sect. [HI It is to be found 
in example_f 90/tabulation_example_streainlined.f90. A program with same interface 
and same output but in F77, is to be found in example_f 77/tabulation_example . f . 

prograun tabulatioii_excLmple_streamlined 
use hoppet_vl 

! ! if using LHAPDF, rename a couple of hoppet functions which 
! ! would otherwise conflict with LHAPDF 

!use hoppet_vl, EvolvePDF_hoppet => EvolvePDF, InitPDF_hoppet => InitPDF 
implicit none 

real(dp) :: dy, ymax, dlnlnQ, Qmin, Qmax, muR_Q 

real (dp) :: asQ, QOalphas, QOpdf 

real(dp) :: mc,mb,mt 

integer : : order, nloop 

! ! holds information about the grid 

type(grid_def ) : : grid, gdarray(4) 

! ! hold results at some x, Q 

real(dp) :: Q, pdf _at_xQ(-6 : 6) 

real(dp), parameter :: heralhc_xvals(9) = & 

& (/le-5_dp , le-4_dp , le-3_dp , le-2_dp , . l_dp , . 3_dp , . 5_dp , . 7_dp , . 9_dp/) 
integer : : ix 

! set up parameters for grid 
order = -6 
ymax = 12.0_dp 
dy = 0. l_dp 

! set up the grid itself — we use 4 nested subgrids 
call InitGridDef (gdarray(4) ,dy/27.0_dp,0.2_dp, order=order) 
call InitGridDef (gdarray(3) ,dy/9.0_dp,0.5_dp, order=order) 
call InitGridDef (gdarray(2) ,dy/3.0_dp, 2. 0_dp, order=order) 
call InitGridDef (gdarray(l) ,dy, ymax , order=order) 

call InitGridDef (grid,gdarray (1 : 4) ,locked= . true . ) 

! Streamlined initialisation 

Qmin=l_dp 

Qmax=28000_dp 

dlnlnQ = dy/4.0_dp 

nloop = 3 

call hoppetStartExtendedCymax , dy , Qmin , Qmax , dlnlnQ , nloop , & 
& order, factscheme_MSbar) 

! Set heavy flavour scheme 
mc = 1.414213563_dp 
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mb = 4 . 5_dp 
mt = 175.0_dp 

call hoppetSetVFN(mc , mb, mt) 

! Set parameters of running coupling 
asQ = 0.35_dp 
QOalphas = sqrt(2.0_dp) 
muR_Q = 1.0_dp 

QOpdf = sqrt(2.0_dp) ! The initial evolution scale 
! Normal evolution 

call hoppetEvolve(asQ , QOalphas, nloop,muR_Q,& 
& LHAsub, QOpdf) 

! Uncomment to perform cached evolution 

! call hoppetPreEvolveCasQ, QOalphas, nloop, muR_Q, QOpdf ) 

! call hoppetCachedEvolve(LHAsub) 

! get the value of the tabulation at some point 
Q = 100.0_dp 
write (6, ' (a) ') 

write(6, ' (a,f8.3,a) ') " Evaluating PDFs at Q = ",Q," GeV" 

write(6, ' (a5,2al2,al4,al0,al2) ' ) "x",& 

& "u-ubar" , "d-dbar" , "2(ubr+dbr) " , "c+cbar" , "gluon" 
do ix = 1, size(heralhc_xvals) 

call hoppetEval(heralhc_xvals(ix) ,Q,pdf _at_xQ) 
write(6, ' (es7. I,5esl2.4) ') heralhc_xvals(ix) , & 
& pdf _at_xQ(2)-pdf _at_xQ(-2) , & 
& pdf _at_xQ(l)-pdf _at_xQ(-l) , & 
& 2*(pdf _at_xQ(-l)+pdf _at_xQ(-2) ) , & 
& (pdf_at_xQ(-4)+pdf_at_xQ(4)) , & 
& pdf_at_xQ(0) 

end do 

contains 

subroutine LHAsub (x,Q, pdf ) 

! Scune as in the previous example progremi 
end subroutine LHAsub 

end progrcun tabulation_example_streamlined 

A. 3 Accessing tables 

As has been mentioned in Sect. 17.21 there are several ways to define, initiahse and access 
a user defined table, for example as in the example below: 

module external_table_module ! common location for your table 
use hoppet_vl 
implicit none 
type (pdf _table) : : table 
type(dglap_holder) : : dh 



48 



type(runiiing_coupling) : : coupling 
type(grid_def ) : : grid ! Optional 
end module external_table_module 

Note that one might need to include in the common module objects like dh or coupling, 
since these are required in other procedures: 

subroutine A 

use external_table_module 

call PreEvolvePdf Table(table, QO, dh, coupling) 
end subroutine A 

subroutine B 

use external_table_module 

call EvolvePdfTable(table,pdf 0) 
end subroutine B 

subroutine C 

use external_table_module 

call EvalPdfTable_xQ(table, x, Q, pdf) 
end subroutine C 

A more detailed description of the above technique can be found in a second example 
program, called tabulation_exainple_2 . f 90, which is available at the example_f90 direc- 
tory. It generates exactly the same output as the previous example program, however it 
provides an example of how to table from external procedures, as explained in 

Sect. Ol 

B HOPPET reference guide 

In this section we present the HOPPET reference guide, a summary of the most impor- 
tant modules in the package with the corresponding description, both the the streamlined 
interface. Table H] and for the general interface. Table [5l 

C Initialisation of grid quantities 

As has been mentioned in Sect. 14.21 there exist several ways of setting a grid quantity. 
In this Appendix we describe the most important methods for initialising a grid quantity, 
which we take to be a parton distribution. 

There are a number of ways of setting a grid quantity. Suppose we have a subroutine 

subroutine exajiiple_gluon(y ,g) 

use types ! ! defines "dp" (double precision) kind 

implicit none 

real (dp), intent (in) :: y 



49 



STREAMLINED INTERFACE 



METHOD 


DESCRIPTION 


Initialisation 




hoppetStart (dy , nloop) 


Sots up a compound grid with spacing in In l/x of dy at small x, 
extending to = 12 and numerical order = — 5. 
The Q range for the tabulation will be 1 GeV < Q < 28 TeV, 
dlnlnQ=dy/4 and the factorisation scheme is MS 


hoppetStartExtendedCymax , dy , Qmin , 
Qmax , dlnlnQ , nloop , order , fact scheme) 


More general initialisation 


hoppetSetFFN(f ixed_nf ) 
hoppetSetVFN(mc , mb, mt) 


Set heavy flavour scheme 


alphas = hoppetAlphaS(Q) 


Accessing the coupling 


Normal evolution 




hoppetEvolve (asQ , QOalphas , 
nloop , muR_Q , LHAsub , QOpdf ) 


PDF evolution: specifies the coupling asQ at a scale QOalphas, 
the number of loops for evol., nloop, 
the ratio (muR_Q) of ren. to fact, scales, 
the name of a subroutine LHAsub with an LHAPDF-like interface 
and the scale QOpdf at which one starts the PDF evolution 
Note: LHAsub only called at scale QOpdf 


hoppetEval(x , Q , f ) 


On return, f (-6:6) contains all flavours of the PDF set 
(multiplied by x) at the given x and Q values 


Cached evolution 




hoppetPreEvolveCasQ , QOalphas , 
nloop, muR_Q, QOpdf) 


Preparation of the cached evolution 


hoppetCachedEvolve (LHAsub) 


Perform cached evolution with the initial condition 
at QOpdf from a routine LHAsub with LHAPDF-like interface 
Notice LHAsub only called at scale QOpdf 


hoppetEval(x , Q , f ) 


On return, f (-6:6) contains all flavours of the PDF set 
(multiplied by x) at the given x and Q values 
[as for normal evolution] 



Table 4: Reference guide for the streamlined interface. Note that hoppet should be 
supplied with and returns parton densities multiplied by x. 
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GENERAL INTERFACE 



TYPES 


DESCRIPTION 


type (grid_def ) :: 


grid 


X— space grid definition 


real (dp), pointer :: 


gluon( : ) 


Holds a 'grid quantity' (e.g. gluon PDF) 


real (dp), pointer :: 


PDFset(: , :) 


Grid representation of a (13-flavour) PDF set 


type (grid_conv) : 


: Pgg 


Convolution operator {i.e. splitting function) 


type (split_mat) : 


: Pmat 


Splitting matrix (with full flavour structure) 


type (mass_threshold_mat) 


: : MTM.NNLO 


Heavy quark mass-threshold matrix 


type (dglap_holder) : 


: dglapji 


DGLAP holder (i.e. all splitting and mass-threshold matrices) 


type (running_coupling) 


: ; coupling 


Running coupling 


type (evln_operator ) 


: : evop 


Evolution operator (linked list of split & mass-threshold matrices) 


type (pdf .table) :: 


table 


PDF set tabulated in x &c Q 



METHOD 


DESCRIPTION 


Initialisation 




InitGr idDef (gr id , dy=0 . l_dp , ymax=10 . 0_dp , order=3) 


Initialise a grid definition 


InitGridDef (grid, subgr ids ( : ) , locked= . true . ) 


Combine subgrids into a single grid 


AllocGridQuant (grid , gluon) 


Allocate memory for a grid quantity (e.g. gluon PDF) 


InitGr idQuant (gr id , gluon, exajnple_gluon_fn) 


Initialisation of a grid quantity (e.g. gluon PDF) 


AllocPDF (grid, pdf set) 


Allocate memory for a 13-fiavour PDF set 


InitPDF.LHAPDF (gr id , pdf set , LHAsub , Q) 


Initialisation of a (13-flavour) PDF set from LHAPDF 
type routine. Note: it only calls LHAsub at the scale Q 


InitGr idConv (grid , Pgg , Pgg_f unc) (*) 


Initialisation of a convolution operator 


InitDglapHolder (grid, dglap_h, factscheme, nloop) (*) 


Initialisation of a dglap_holder type 


lnitRunningCoupling(coupling [, alfas] [, Q] 
[, nloop] [, fixnf] [, quarkjnasses] ) (*) 


Initialisation of a running_coupling type 


AllocPdf Table (grid, table, Qmin, Qmax 
[, dlnlnQ ] [, InlnQ.order ] [, f reeze_at_Qmin] ) (*) 


Allocate space for a pdf .table type 


Evaluation & manipulation 




EvalGr idQuant (gr id , gluon , y ) 


Evaluation of a grid quantity at y = Inl/x 


Pgg. conv. gluon 


Convolution of a splitting function with a (1-flav) PDF 


Pmat . conv . PDFset 


Convolution of a splitting matrix with a PDF set 


SetToConvolution(Pab , Pac , Pcb) 


Convolution of splitting functions, Pab = Pac Pcb 


EvalPdfTable.yQ (table, y, Q, pdf) 


Evaluate the 13 flavours, pdf (-6:6), of a tabulated 
PDF set at y and Q 


Evolution 




CopyHumanPdf ToEvln(nf _lcl , pdfjiuman, pdf_evln) 


Transform PDF set from human to evln representation 


GetPdf Rep (pdf set) 


Check PDF set representation 


EvolvePDF(dglap_h, initial_pdf set , coupling, 
Q_init,Q_end , [, muR.Q] [, nloop] [,du]) 


Evolution of an initial condition for a PDF set 
at Q_init to Q_end 


EvolvePdf Table (table, QO, initial.pdf set , dglap_h, 
coupling [, muR_Q] [, nloop] [, untie_nf]) 


Fill a table starting from an initial condition, 
initial_pdf set, at the scale QO 



Table 5: Reference guide for the general interface. The upper table describes the main 
derived types defined in hoppet. The lower table summarises some of the most relevant 
methods. Note that arguments between [ . . . ] are optional. For grid quantities and PDF 
sets the user must explicitly make memory allocation calls; in other cases (marked with 
a (*)), initialisation routines automatically allocate the memory. For all types, allocated 
memory may be freed with a call to the Delete ( . . . ) subroutine. 
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real (dp), intent (out) :: g 
real (dp) : : x 



X = exp(-y) 

g = 1.7_dp * x**(-0.1_dp) * (l-x)**5 
end subroutine example_gluon 

Then we can call 

call lnitGridQuantSub(grid,gluon,excunple_gluon) 

to initialise gluon with a representation of the return value from the subroutine exainple_gluon. 

An alternative way is to make use of functions xValues or yValues that respectively 
return the x or y values of all points on the grid: 

real(dp), pointer :: gluon, xvals 
call AllocGridQuant(grid, gluon) 
call AllocGridQuant(grid, xvals) 
xvals = xValues(grid) 

gluon = 1.7_dp * xvals** (-0 . l_dp) * (l-xvals)**5 
deallocate (xvals) 

Though more laborious insofar as one has to worry about some extra allocation and deal- 
location, it has the advantage that one no longer has to write a separate subroutine. 

Finally, there is an option to initialise a multi-flavour PDF grid with a subroutine with 
the same format as evolvePDF from the LHAPDF library. This option works as follows: 

real (dp) , pointer :: pdf_set(:,:) 

real (dp) : : y,Q 

real(dp) :: pdf _at_y(-6 : 6) 

Q=2 ! Initial scale 

! Initialise with LHAPDF-like routine 
call AllocPDF(grid,pdf _set) 

call InitGridQuantLHAPDF(grid, pdf_set, LHAsub, Q) 

! Evaluate the multi-flavor pdf at y 

pdf_at_y = EvalGridQuant(grid,pdf _set ( : , -6 : 6) ,y) 

where an example of the LHAsub subroutine can be found in Sect. 15.1.11 

D NNLO splitting functions 

A remark concerning NNLO splitting functions: the exact NNLO splitting functions de- 
rived by Moch, Vermaseren and Vogt [2T], [22] involve long (multi-page) expressions in terms 
of harmonic polylogarithms of up to weight 4. Very conveniently, refs. [211 122] provide the 
expressions directly in terms of Fortran code. The harmonic polylogarithms can be evalu- 
ated using the hplog package of Gehrmann and Remiddi [13], a copy of which is included 
with the HOPPET package. 
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The initial integrations needed to create the splitjnat objects for the exact NNLO 
sphtting functions for the full range of take of the order of minutes. Since currently 
there is no option of storing the splitting matrices in a file, this can be a bit bothersome. 
So instead, by default, the program uses the approximate, parametrised NNLO splitting 
functions also provided in [211 [22]. The parametrised splitting functions are guaranteed to 
be accurate to within 0.1% — in practice since they come in relatively suppressed by two 
powers of the impact on the evolution tends to be of the order of a 10~^ relative effect 

The user can choose whether to obtain the exact or parametrised NNLO splitting 
functions using the following calls (to be made before initialising the splitting matrices) 

integer : : splitting_varicLnt 

call dglap_Set_niilo_splitting(splitting_variant) 

with the following variants defined (as integer parameters) in the module dglap_choices: 
nnlo_splitting_exact 

nnlo_splitting_param [default] 
nnlo_splitting_Nf itav 
nnlo_splitting_Nf iterrl 
nnlo_splitting_Nf iterr2 

The last 3 are the parametrisations based on fits to reduced moment information carried 
out in [301 EI]. Though at the time they represented a valuable (and much used) step on 
the way to full NNLO results, nowadays their interest is mainly historical. 

Note that only for the nnlo_splitting_exact can the colour constants be varied (with 
the caveat about d°'^'^dabc , as is described in Sect. 15.3.11) . For the other options the NNLO 
splitting functions have been computed using the QCD values for the colour factors. 

E Useful tips on Fortran 95 

As Fortran 95's use in high-energy physics is not as widespread as that of other languages 
such as Fortran 77 and C-I--I-, it is useful to summarise some key novelties compared to 
Fortran 77, as well as some points that might otherwise cause confusion. For further 
information the reader is referred both to books about the language such as [H] and to 
web resources [15] . 

Free form. Most of the code in the hoppet package is in free-form. The standard exten- 
sion for free- form form files is . f 90. There is no requirement to leave 6 blank spaces before 
every line and lines can consist of up to 132 characters. The other main difference relative 
to f77 fixed form is that to continue a line one must append an append an ampersand, &, to 
the line to be continued. One may optionally include an ampersand as the first non-space 
character of the continuation line. 

For readability, many of the subprogram names in this documentation are written with 
capitals at the start of each word. Note however that free- form Fortran 95, like its fixed- 
form predecessors, is case insensitive. 
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Modules, and features relating to arrays. Fortran 95 allows one to package variables 
and subroutines into modules 

module test_module 
implicit none 
integer : : some_integer 
contains 

subroutine print_array(array) 

integer, intent(in) :: array(:) ! size is known, first element is 1 

! intent (in) == array will not be changed 

integer : : i , n 

n = size (array) 
do i = 1 , n 

print *, i, array (i) 
end do 

end subroutine hello_world 
end module test_module 

The variable some_integer and the subroutine print_array are invisible to other routines 
unless they explicitly use the module as in the following example: 

progremi test_program 
use test_module 
implicit none 

integer :: arrayl(5), array2(-2:2) 
integer : : i 



some_integer = 5 
array 1 = 

array2(-2:0) = 9£ 



set the variable in test_module 

set all elements of arrayl to zero 

set elements 1..3 of array2 to equal to 3. 



array2(l:2) = 2*array2(-l : 0) ! elements -2..0 equal twice elements -1..0 

print *, "Printing array 1" 
call print_array(arrayl) 
print *, "Printing array 2" 
call print_array(array2) 
end progreim test_prograin 

Constants can be assigned to arrays (arrayl) or array subsections (array2(-2 : 0)), arrays 
can be assigned to arrays of the same size (as is done for array2 (-2 : 0) ) and mathematical 
operations apply to each element of the array (as with the multiplication by 2). 

When arrays are passed to function or subroutine that is defined in a used module, 
information about the size of the array is passed along with the array itself. Note however 
that information about the lower bound is not passed, so that for both arrayl and array2, 
print_array will see arrays whose valid indices will run from 1 ... 5. Thus the output from 
the program will be 

Printing array 1 

1 

2 

3 

4 
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5 

Printing array 2 

1 99 

2 99 

3 99 

4 198 

5 198 

If print_array wants array to have a different lower bound it must specify it in tlie 
declaration, for example 

integer, intent(in) :: array(-2:) ! size is known, first element is -2 

While it may initially seem bizarre, there are good reasons for such behaviour (for example 
in allowing a subroutine to manipulate multiple arrays of the same size without having to 
worry about whether they all have the same lower bounds). 

Dynamic memory allocation, pointers. One of the major additions of f95 compared 
to f77 is that of dynamic memory allocation, for example with pointers 

integer, pointer :: dynamic_array( : ) 

allocate (dynamic_array (-6 : 6) ) 

! . . work with it . . 

de al 1 o c at e ( dynami c _ array) 

This is fundamental to our ability to decide parameters of the PDF grid(s) at run-time. 
Pointers can be passed as arguments to subprograms. If the subprogram does not specify 
the pointer attribute for the dummy argument 

subroutine xyz(dummy_array) 

integer, intent(in) :: duininy_array( : ) 

then everything behaves as if the argument were a normal array (e.g. the default lower 
bound is 1). Alternatively the subroutine can specify that it expects a pointer argument 

subroutine xyz (dummy_pointer_array) 

integer, pointer :: dummy _pointer_array( : ) 

In this case the subroutine has the freedom to allocate and deallocate the array. Note also 
that because a pointer to the full array information is being passed, the lower bound of 
dummy _pointer_arr ay is now the same as in the calling routine. Though this sounds like a 
technicality, it is important because a corollary it that a subroutine can allocate a dummy 
pointer array with bounds that are passed back to the calling subroutine (we need this for 
the flavour dimension of PDFs, whose lower bound is most naturally —6). 

Note that in contrast to C/C++ pointers, F95 pointers do not explicitly need to be 
dereferenced — in this respect they are more like C++ references. To associate a pointer 
with an object, one uses the => syntax: 

integer, target : : target_object(10) 
integer, pointer :: pointer_to_object( : ) 
pointer_to_object => target_object 

pointer_to_object(l : 10) = ! sets target_object(l : 10) 
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One notes that the object that was pointed to had the target attribute — this is mandatory 
(unless the object is itself a pointer). 

Derived types. Another feature of F95 that has been heavily used is that of derived 
types (analogous to C's struct): 

type pair 

integer first, second 
end type pair 

Variables of this type can then be created and used as follows 

type (pair) : : pair_object, £m.other_pair_object 

pair_object%f irst = 1 

pair_object%second = 2 

another _pair_object = pair_object 

print *, another_pair_objecty„second 

where one sees that the entirety of the object can be copied with the assignment (=) 
operator. Note that many of the derived types used in hoppet contain pointers and when 
such a derived type object is copied, the copy's pointer just points to the same memory 
as the original object's pointer. This is sometimes what you want, but on other occasions 
will give unexpected behaviour: for example splitting function types are derived types 
containing pointers, so when you assign one splitting function object to another, they end 
up referring to the same memory, so if you multiply one of them by a constant, the other 
one will also be modified. 

Operator overloading While assignment behaves more or less as expected by default 
with derived types (it can actually be modified if one wants to), other operators do not 
have default definitions. So if one wants to define, say, a multiplication of objects one may 
associate a function with a given operator, using an interface block: 

module test_module 

interface operator(*) ! provide access to dot_pairs through 

module procedure dot_pairs ! the normal multiplication symbol 
end interface 

interface operator (. dot . ) ! provide access to dot_pairs through 

module procedure dot_pairs ! a specially named operator 
end interface 
contains 

integer function dot_pairs(pairl , pair2) 
type (pair), intent (in) :: pairl, pair2 

dot_pairs = pairlZf irst*pair2yof irst + pairl7.second*pair2%second 
end function dot_pairs 
end module 

given which we can then write 

integer : : i 

type (pair) :: pairl, pair2 

[. . . some code to set up pair values . . .] 
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! now multiply them 
i = pairl * pair2 

i = pairl .dot. pair2 ! equivalent to previous statement 

Since the the multiphcation operator (*) aheady exists for all the default types, by defin- 
ing it for a new type we have overloaded it. Note that there are some subtleties with 
precedences of user-defined operators: operators (like *) that already exist have the same 
precedence as they have is usual operators; operators that do not exist by default (.dot) 
have the lowest possible preference, so, given the above definitions, 

i = 2 + pairl * pair2 ! legal 

i = 2 + pairl .dot. pair2 ! illegal, means: (2+pairl) .dot .pair2 
i = 2 + (pairl .dot. pair2) ! legal 

where the second line is illegal because we have not defined any operator for adding an 
integer and a pair. Similarly care is needed when using the hoppet's operator .conv.. 

Floating point precision: A final point concerns floating point variable types. Through- 
out we have used definitions such as 

real (dp), pointer :: pdf ( : , : ) 

and written numbers with a trailing _dp 

param = 1 . 7_dp 

Here dp is an integer parameter (defined in the types module and accessible also through 
the hoppet_vl module), which specifies the kind of real that we want to define, specifically 
double precision. We could also have written double precision everywhere, but this is 
less compact, and the use of a kind parameter has the advantage that we can just modify its 
definition in one point in the program and the precision will be modified everywhere. (Well, 
almost, since some special functions are written in Fortran 77 using double precision 
declarations and do their numerics based on the assumption that that truly is the type 
they're dealing with). 

Optional and keyword arguments A feature of F95 that helps simplify user interfaces 
is that of optional and keyword arguments. Suppose we have 

subroutine hello(name, prefix, count) 

character(len=*) , intent(in) :: name, prefix 

integer, optional, intent(in) : : count 
end subroutine hello 

Here the count argument is optional meaning that it need not be supplied — if it is absent 
the subroutine is should behave sensibly all the same. Thus one can call the subroutine as 

call hello(name, prefix) 

call hello (name, prefix, count) 

Keyword arguments are useful if one doesn't want to remember the exact order of a long 
list of arguments (or if one wants to specify just one of several optional arguments). For 
example 
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call hello (rLajne=iiame, pref ix=pref ix) 
call helloCpref ix=pref ix, naine=nEiine) 

will do the same thing. 
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